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We investigate properties of material ejected dynamically in the merger of black hole-neutron 
star binaries by numerical-relativity simulations. We systematically study the dependence of ejecta 
properties on the mass ratio of the binary, spin of the black hole, and equation of state of the neutron- 
star matter. Dynamical mass ejection is driven primarily by tidal torque, and the ejecta is much 
more anisotropic than that from binary neutron star mergers. In particular, the dynamical ejecta is 
concentrated around the orbital plane with a half opening angle of 10°-20° and often sweeps out only 
a half of the plane. The ejecta mass can be as large as ~ O.IM©, and the velocity is subrelativistic 
with ~ 0.2-0.3c for typical cases. The ratio of the ejecta mass to the bound mass (disk and fallback 
components) is larger, and the ejecta velocity is larger, for larger values of the binary mass ratio, 
i.e., for larger values of the black-hole mass. The remnant black hole-disk system receives a kick 
velocity of 0(100) kms“^ due to the ejecta linear momentum, and this easily dominates the kick 
velocity due to gravitational radiation. Structures of postmerger material, velocity distribution of 
the dynamical ejecta, fallback rates, and gravitational waves are also investigated. We also discuss 
the effect of ejecta anisotropy on electromagnetic counterparts, specifically a macronova/kilonova 
and synchrotron radio emission, developing analytic models. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 


I. INTRODUCTION 

Coalescences of black hole-neutron star binaries are 
one of the most promising gravitational-wave sources for 
ground-based laser-interferometric detectors along 
with those of binary neutron stars and binary black holes. 
The sensitivity of these detectors will reach a sufficiently 
high level in the coming years to detect gravitational 
waves from compact binary coalescences more often than 
once a year mg- The first direct detection of gravita¬ 
tional waves must have a dramatic impact on fundamen¬ 
tal physics. Furthermore, gravitational waves from bi¬ 
naries involving neutron stars will tell us neutron-star 
properties like the radius, compactness, and tidal de- 
formability. Knowledge of neutron-star properties will 
allow us to constrain the equation of state of nuclear- and 
supranuclear-density matter, and therefore gravitational 
waves will also give us valuable information on nuclear 
physics. 

Simultaneous detection of electromagnetic radiation 
from compact binary mergers, i.e., electromagnetic coun¬ 
terparts to gravitational waves, is eagerly desired liH]. 
It will support gravitational-wave detection and en¬ 
hance scientific returns from each coalescence event. 
For example, source localization on the celestial sphere 
is much more accurate with electromagnetic instru¬ 


ments than with gravitational-wave detector networks 
[5] . Gravitational-wave data analysis benefits from accu¬ 
rate localization by solving degeneracy between the sky 
location and other amplitude parameters such as the lu¬ 
minosity distance. Accurate localization of the source is 
also indispensable to find its host galaxy and to deter¬ 
mine the cosmological redshift. By combining this in¬ 
formation, the luminosity distance-redshift relation will 
be derived without relying on the cosmic distance lad¬ 
der}^ and we will obtain a novel method to test cosmo¬ 
logical models [10]. Besides, the effective sensitivity of 
a gravitational-wave detector would be improved if we 
could know the coalescence time and/or sky location of 
a binary from electromagnetic counterparts HH. 

Among the candidates of electromagnetic counter¬ 
parts, a short-hard gamma-ray burst and its afterglow 
are vigorously studied both theoretically and observa- 
tionally (see Refs. nana for reviews). While whether 
compact binary coalescences can really drive short-hard 
gamma-ray bursts is still an open question, future simul¬ 
taneous detection with gravitational-wave chirp signals 
will prove this hypothesis. Prompt emission is so bright 


^ See also Ref. for an alternative approach free from electro¬ 
magnetic observation. 
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that it can be easily detected by gamma-ray satellites 
within the horizon distance of gravitational-wave detec¬ 
tors. Accurate localization is possible if an associated 
afterglow is observed at longer wavelengths. Short-hard 
gamma-ray bursts do not, however, always serve as coun¬ 
terparts to gravitational waves because of their presum¬ 
ably jetlike geometry. If the typical jet opening angle 
is < 10° as suggested by jet-break observations [Til ITS] , 
the fraction of gravitational-wave events accompanied by 
observable short-hard gamma-ray bursts will be a few 
percent at best. 

In recent years, electromagnetic counterparts have 
been getting a lot more attention, and many isotropic 
emission models are studied. Most of the proposed mod¬ 
els require the ejection of unbound material from the 
binarym where examples include a macronova/kilonova 
powered by decay heat of unstable r-process elements 
PD1 - E5] and nonthermal radiation from electrons acceler¬ 
ated at blast waves between the ejecta and interstellar 
medium [241 [25]. Possible emission from the ejecta will 
be isotropic if the ejecta has spherical geometry and/or a 
subrelativistic velocity. Such a “47r-counterpart” is ideal 
for followup observations, because it will accompany a 
majority of gravitational-wave events unlike beamed ra¬ 
diation of short-hard gamma-ray bursts. 

Despite its 40-years-old history in theoretical astro¬ 
physics HSIIIT], mass ejection from the compact binary 
merger is a young research topic in numerical relativ¬ 
ity. Most of the previous black hole-neutron star binary 
simulations in numerical relativity were performed aim¬ 
ing at deriving gravitational waves in the late inspiral 
and merger phases and at clarifying properties of rem¬ 
nant accretion disks formed after the tidal disruption of 
neutron stars (see Ref. [55] and references therein for ear¬ 
lier work). Mass ejection has not been studied in detail 
compared to these topics in full general relativity [5i- 
[83] . whereas a substantial effort to clarify mass ejection 
has been made in simulations performed in Newtonian 
gravity or approximate general relativity [54H55] (see also 
Refs. [57M5] 1. It is pointed out that dynamical ejecta 
from binary neutron star mergers become less massive 
and more isotropic in full general relativity [44] or the 
conformal flatness approximation |45| than in Newtonian 
gravity [35l [33 [47] . The difference due to realism of the 
gravity should be most pronounced when a black hole is 
involved as already suggested by existing work. Thus, it 
is natural that numerical relativity is vital to study mass 
ejection from the black hole-neutron star binary merger. 

In this study, we perform simulations of black hole- 
neutron star binary mergers using numerical-relativity 
code SACRA [35] and investigate dynamical mass ejec¬ 
tion extending our preceding work |49| . In particular, 
we focus on kinematical properties of dynamical ejecta 


such as the mass and velocity. Compared to our pre¬ 
vious simulations [29] , we adopt large computational do¬ 
mains to track long-term evolution of the ejecta. Because 
the dynamical ejecta has a velocity comparable (a few 
tens of percent) to the speed of light as shown in this 
study, the large computational domains are essential for 
the reliable estimation of ejecta properties. We also im¬ 
prove the treatment of artificial atmosphere (inevitable 
in conservative hydrodynamic schemes) from our previ¬ 
ous work [29l 1351 [5nU53] and confirm that characteristic 
quantities of dynamical ejecta such as the mass and ve¬ 
locity depend only weakly on the atmosphere. We do 
not, however, study disk winds expected to be driven by 
unincorporated physics. 

This paper is organized as follows. Section [ll] describes 
our models of black hole-neutron star binaries including 
neutron-star equations of state. Our numerical meth¬ 
ods are also described, and diagnostics of simulations are 
presented with particular emphasis on the ejecta defined 
as unbound material. Numerical results of the simula¬ 
tions are presented in Sec. El After briefly reviewing 
the merger dynamics in Sec. Ill A[ dynamical mass ejec¬ 
tion processes are described in Sec. jlllB] The depen¬ 
dence of characteristic quantities on binary parameters 
is discussed in Sec. Ill C and the material structure is 
investigated in Sec. HID We also study fallback mate¬ 


rial, remnant black hole-disk systems, and gravitational 
waves in Secs. jlll E[ [IITfI and [III G] respectively. Possi¬ 
ble electromagnetic counterparts from black hole-neutron 
star binaries are discussed based on the results of simu¬ 
lations in Sec. HYl Specifically, Sec. jlV Aj describes the 
macronova/kilonova and Sec. IV B[ describes synchrotron 
radio emission from accelerated electrons. Section E] is 
devoted to a summary. Numerical values of characteristic 
ejecta quantities derived by simulations are summarized 
in Table III Readers interested primarily in electromag¬ 
netic counterparts should read Sec. |IV[ of which impor¬ 
tant results are described in Ref. [45] . 

Notational conventions are summarized as follows. 
Throughout this paper, we adopt geometrical units in 
which G = c = 1, where G and c are the gravitational 
constant and speed of light, respectively. Exceptionally, 
c is sometimes inserted for clarity when we discuss the 
velocity of the ejecta or fluid element. Greek and Latin 
indices denote the spacetime and space components, re¬ 
spectively. The black-hole mass, neutron-star gravita¬ 
tional mass, and neutron-star circumferential radius in 
isolation are denoted by Mbh, AIns, and i?NSj respec¬ 
tively. The dimensionless spin parameter of the black 
holej^ total mass of the system at an infinite separa¬ 
tion, mass ratio, and compactness of the neutron star 
are defined as y = mo = Mbh + Mns, 

Q = Mbh/Mns, and C = Mns/Rns, respectively, where 


Precursor emission may not require mass ejection [16ttl9| . and 
we do not consider them in this study. 


^ In our previous work |29M52lf5^ . this parameter is denoted as a. 
We change the convention, because a is sometimes reserved for 
the specific spin angular momentum, 5 'bh/-^bh = X-^BH- 
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S'bh is the black-hole spin angular momentum. 

II. NUMERICAL METHOD 
A. Zero-temperature equation of state 

We model equations of state of zero-temperature 
neutron-star matter by piecewise polytropes [M]. Neu¬ 
tron stars in compact binaries right before the merger will 
be cold enough to be modeled by zero-temperature equa¬ 
tions of state (see, e.g., Ref. [55]). However, the realistic 
equation of state of neutron-star matter is not known 
precisely yet. Therefore, it is necessary to adopt vari¬ 
ous equations of state systematically to span a plausible 
range of neutron-star properties. Piecewise polytropes 
are suitable for this purpose, because those with one and 
three pieces for crust and core regions, respectively, are 
known to be able to approximate nuclear-theory-based 
equations of state accurately with a small number of pa¬ 
rameters |S3]. Following Ref. |^, we employ piecewise 
polytropes of the form 

P{p) = {p,< p< Pi+i), (1) 

where P is the pressure and p is the rest-mass density, 
with i G {0,1,2,3} in this study. It is always assumed 
that po = 0 and p 4 —>■ oo. We fix parameters for the 
lowest-density, crust region to be 

Fo = 1.356 923 95, (2) 

Ato = 3.998 736 92 X 10"® (gcm-^)^”^”. (3) 

We further set p 2 = 10^"^'’^ gcm“® « 5.0 x 10^^gcm“® 
and ps = 10 ^® ° gcm“® to reduce the number of free pa¬ 
rameters. Requiring the continuity of P{p), each piece- 
wise polytrope is characterized by four parameters. We 
choose the free parameters to be the pressure at p 2 , de¬ 
noted by P 2 = P{p 2 ), and adiabatic indices for the core, 
{ri,r2,r3}. 

Table |l] lists parameters of piecewise polytropes 
adopted in this study as well as neutron-star properties 
computed using them. The naming convention and pa¬ 
rameters follow Ref. |S1|. APR4 is computed by a 
variational method incorporating three-nucleon interac¬ 
tions and relativistic boost corrections. This equation of 
state gives the smallest radius of a 1 . 35 M 0 neutron star, 
R1.35 = 11.1km, and thus APR4 is the softest equation 
of state among those adopted in this study. Accordingly, 
tidal disruption is less pronounced for neutron stars mod¬ 
eled with APR4 than those modeled with the other equa¬ 
tions of state. ALF2 m is a hybrid equation of state 
obtained combining a nucleonic, APR-type equation of 
state at low density and a quark-matter equation of state 
with quantum chromodynamics corrections at high den¬ 
sity. H4 [58l |59j is computed by relativistic mean-field 
theory incorporating hyperons with the stiffest possible 
parameters (at the time). MSI [50] is also derived by rel¬ 
ativistic mean-field theory for nucleonic matter and gives 


R1.35 = 14.4 km, which is the largest value in this study. 
Thus, MSI is an extreme example with which tidal dis¬ 
ruption occurs most violently. 

In practice, a very-high-density regime is not rele¬ 
vant to black hole-neutron star binary coalescences as 
far as canonical-mass neutron stars with M^g « 1.35 Mq 
EH ED are concerned. The reason for this is that the 
maximum rest-mass density of the system, i.e., the cen¬ 
tral density of the neutron star and maximum density in 
the remnant accretion disk, is a decreasing function of 
time except for subdominant oscillations. The rest-mass 
density at the center of an isolated 1.35Mq neutron star, 
P 1 . 35 , never exceeds p^ for the equations of state adopted 
in this study (see Table |T]), and hence Fa never plays a 
role in black hole-neutron star binary coalescences. For 
MSI, even r 2 is irrelevant, because pi .35 is lower than 
P 2 0 This situation is in stark contrast to that of binary 
neutron star coalescences, which depend crucially on the 
high-density regime of the equations of state. 

In this study, we regard quantities associated with 
1.35Mq neutron stars as characteristic quantities of the 
equation of state rather than the maximum mass Mmax, 
which is sensitive to the behavior of matter at high den¬ 
sity. Table|l]shows the baryon rest mass M*, compactness 
C, quadrupolar tidal Love number k |63LI64) . and dimen¬ 
sionless quadrupolar tidal deformability A = (2/3)/cC“® 
of a 1.35Mq neutron star, in addition to i?i, 35 , ^ 1 , 35 , and 
Almax- Note that all the equations of state can support 
~ 2 M 0 neutron stars and satisfy constraints from ob¬ 
servations of massive pulsars |65l [ 66 | , and thus they are 
possible candidates of the realistic equation of state. 


B. Initial condition 

We adopt quasiequilibrium states of black hole-neutron 
star binaries as initial data of our simulations in the same 
manner as Refs. ED ED ESj • Here, we briefly describe the 
computational method of quasiequilibrium states, and 
the details are found in Refs. EDETj. Numerical com¬ 
putations are performed using the multidomain spectral 
method library LORENE [55] . 

We solve a subset of the Einstein equation and the 
hydrostatic equilibrium equations assuming the exis¬ 
tence of helical symmetry. Hamiltonian and momen¬ 
tum constraints are solved by a mixture of the con¬ 
formal transverse-traceless decomposition [69] and ex¬ 
tended conformal thin-sandwich formulation (^ED im¬ 
posing the spatial conformal flatness, maximal slicing, 
and preservation of them in time. The singularity as¬ 
sociated with the black hole is handled in the punctnre 


This means that two-piecewise polytropes adopted in Refs. |29[ 
1521153| can fully replace four-piecewise polytropes adopted here 
for modeling such a stiff equation of state in simulations of black 
hole-neutron star binary coalescences. 
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TABLE I. Parameters and key ingredients of the adopted equations of state, where P 2 (the pressure at p = 10^'^'^ gcm~^) is 
shown in units of dynecm~^. M^ax is the maximum gravitational mass of a spherical neutron star for a given equation of 
state. i?i.35, P1.35, M*,i. 35 , Cl, 35, fei,35, and A1.35 are the circumferential radius, central rest-mass density, baryon rest mass, 
compactness, quadrupolar tidal Love number, and dimensionless quadrupolar tidal deformability of a 1.35Mq neutron star, 
respectively. 


Model 

logio(P2) 

Li T 2 

r3 

Mniax[MQ] 

R 1.35 (km) 

P 1.35 (gcm-^) 

M*,i,35 [Mq] 

Cl. 35 

fci .35 

A 1.35 

APR4 

34.269 

2.830 3.445 

3.348 

2.20 

11.1 

8.9 X 10"^ 

1.50 

0.180 

0.0908 

323 

ALF2 

34.616 

4.070 2.411 

1.890 

1.99 

12.4 

6.4 X 10^"^ 

1.49 

0.161 

0.120 

734 

H4 

34.669 

2.909 2.246 

2.144 

2.03 

13.6 

5.5 X 10^^ 

1.47 

0.147 

0.115 

1110 

MSI 

34.858 

3.224 3.033 

1.325 

2.77 

14.4 

4.2 X 10^"^ 

1.46 

0.138 

0.132 

1740 


framework m, and thus we obtain initial data of the in¬ 
duced metric jij and extrinsic curvature Kij everywhere 
on the initial hypersurface (except for the exact location 
of the puncture, with which simulation grids are chosen 
not to coincide). The neutron-star matter is modeled by 
a perfect fluid expressed by an energy-momentum tensor 
of the form 

Tpii/ = phu^Uy -j- Pg^i/^ (4) 

where h = 1 + e + {P/p) is the specific enthalpy, e is the 
specific internal energy, and is the fluid four-velocity. 
We further assume that the fluid is in a zero-temperature 
and irrotational state during the computation of the ini¬ 
tial data, and hydrostatic equilibrium configurations are 
obtained by solving the continuity equation and local 
energy-momentum conservation equation |78I - I76] . 

Parameters characterizing a black hole-neutron star 
binary are specified in initial data computations (see 
Refs. [29l|67j for the details). For simplicity, we always 
choose Mns to be a typical mass of observed binary neu¬ 
tron stars, Mns = l.SSM©, in this study [ST1I21]. With 
this choice, the black-hole mass, Mbh, is uniquely de¬ 
termined by the mass ratio, Q, which we regard as an 
independent parameter instead of Mbh itself. We only 
consider cases in which the spin angular momentum of 
the black hole is zero or alimed with the orbital angu¬ 
lar momentum of the binary}^ and thus the spin is fully 
characterized by its dimensionless magnitude, y. The or¬ 
bital angular velocity of a binary LI is determined by a 
force-balance condition at the center of the neutron star 
for a given orbital separation. We use a dimensionless 
orbital angular velocity mofl to characterize the initial 
data rather than the orbital separation. 

C. Dynamical simulation 


evolve the conformal-factor variable LF, conformal metric 
7 y-, conformal connection function T*, extrinsic curvature 
trace K, and conformally weighted traceless part of the 
extrinsic curvature Aij defined by 

W = , f* = (5) 

K = , A,, = 7 - 1/3 , ( 6 ) 

where Cartesian coordinates are adopted. The lapse 
function a and shift vector /3* are evolved by a moving 
puncture gauge condition [SOI [SI] of the form 

(dt - P^dj) a = -2aK, (7) 

{dt - /3* = \b\ (8) 

{dt - R* = {dt - P^d,) r - (9) 

where R* is an auxiliary vectorial variable and rjs is a free 
parameter. Initial data of the lapse function are given by 
a = IT, and the shift vector is initialized as /?* = 0 with 
R® = 0. We adopt ps ~ I/wq in this study. 

Hydrodynamic evolution equations are solved by a 
high-resolution shock-capturing central scheme [55] with 
third-order piecewise parabolic reconstruction (53). We 
evolve the conserved rest-mass density p*, conserved mo¬ 
mentum density p^Ui, and conserved energy density p*e 
defined via 

P 

p* = pa^/pv/ , Ui = hui , e = hau* - 7. (10) 

paw* 

Equations of state adopted in dynamical simulations 
comprise cold and thermal parts. The former is taken 
to be piecewise polytropes described in Sec. m and 
the latter is given by an ideal-gas-like form 


Our numerical simulations are performed using an 
adaptive-mesh-refinement code, SACRA [JS]. The Ein¬ 
stein evolution equations are solved in a Baumgarte- 
Shapiro-Shibata-Nakamura formulation [7i IZnj. We 


^ We will report results of cases in which the black-hole spin an¬ 
gular momentum is inclined with respect to the orbital angular 
momentum in a subsequent paper Ez). 


Pth — (Tth l)p^thj 


( 11 ) 


where the thermal-part specific internal energy is defined 
by eth(p, e) = e-ecoid(p) with ecoid(p) the cold-part spe¬ 
cific internal energy computed by piecewise polytropes. 
Total pressure is given by P = Pcoid(p) + Pth(p, £), where 
Pcoid(p) is computed by piecewise polytropes. We choose 
a fiducial value of Tth to be 1.8 following Ref. [44] (see 
also Ref. [S3]) and als o ad opt 1.6 and 2.0 for selected 
models (see Appendix A3). Note that these values are 
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larger than that adopted in our previous work [2111111 [S3], 
in which Fth is always chosen to be Fq [see Eq. ([^]. 

In our simulations, all the postmerger material is gov¬ 
erned effectively by the same sub-nuclear-density equa¬ 
tion of state irrespective of the adopted piecewise poly¬ 
trope. Specifically, when the rest-mass density falls below 
Pi, the equation of state is given by the sum of the crust 
polytrope and thermal correction, 

P{p < pi, e) = Kop^“ -I- (Fth - l)p£th- (12) 


The values of pi are computed as (kqP? 
for each piecewise polytrope (see Sec. II A) and take 0.9- 
2x lO^^gcm”^. The rest-mass density never exceeds 
these values after tidal disruption of neutron stars. 

An artificial atmosphere has to be set carefully to study 
mass ejection accurately. According to Ref. [H], we put 
an atmospheric density floor of the form 


P*,at — /atP*,o min 



(13) 


where p*_o is tbe maximum (conserved) rest-mass density 
of the initial configuration (see Ref. [15] for our previous 
treatment). We typically choose /at = 10“^^ and rt at = 3 
and vary them for selected models (see Appendix |A 2[ ) . 
The critical radius i?crit is chosen to be T/16, where L 
is the size of the computational domain on one side (see 
below). The atmospheric velocity is set to be zero, and 
the atmospheric pressure is given by zero-temperature 
equations of state. 

The grid structure of SACRA is summarized as follows. 
Computational domains are composed of nested equidis¬ 
tant Cartesian grids, and each grid has (2iV -|- 1,2A^ -|- 
1,A^ -I- 1) points in (x,y,z) directions. The equatorial 
symmetry is imposed on the z = 0 plane. We adopt 
fV = 60 as a fiducial value, with which the neutron-star 
radius is covered by > 50 points in the finest grid. We 
also perform simulations with = 40 and 48 for se¬ 
lected models to c heck the convergence of ejecta proper¬ 
ties (see Appendix |A 1[ ). The outer boundary is a cuboid 
covering {x,y,z) € [—L : L] x [—L : L] x [0 : L], and 
outgoing-wave boundary conditions are imposed except 
for the z = 0 plane. As for the adaptive-mesh-refinement 
grid structure, we prepare Ic coarser nonmoving grids 
and If finer moving grids. Namely, we have + 21 f 
computational grids spanning I,, + If refinement levels, 
which we always choose to be 9 in this study. The non¬ 
moving grids are fixed around an approximate center of 
mass throughout the simulation. One set of the mov¬ 
ing grids follows the black hole, and the other set fol¬ 
lows the neutron star. Starting from the coarsest level as 
I = 0, the Rh level has a grid spacing Axi = L/{2’‘N), 
and we specifically denote the grid spacing at the finest 
level by Ax = L/( 2 *‘=+b-i^/ Finally, time steps of 
all the moving grids are chosen by setting the Courant- 
Friedrichs-Lewy factor to be 0.5, and those of the non¬ 
moving grids are chosen to agree with that of the Icth 
level (i.e., the coarsest moving grid). In other words, the 


Courant-Friedrichs-Lewy factor is given by 0.5/2*= * in 
the nonmoving grids. 


D. Binary model and grid setting 

Table [IT] lists black hole-neutron star binary models 
considered in this study. We name each model after 
the equation of state, mass ratio, and black-hole spin. 
For example, APR4-Q3a75 is a binary modeled with the 
APR4 equation of state, Q = 3, and y = 0.75. Recall 
that A/ns = I. 35 M 0 for all the models. Table |ll] also 
presents the dimensionless initial orbital angular veloc¬ 
ity mo^lo, Arnowitt-Deser-Misner mass Mq, and orbital 
angular momentum of the system Jg. Here, Jg is de¬ 
fined from an Arnowitt-Deser-Misner type integral by 
subtracting the spin angular momentum associated with 
the puncture. 

We take the mass ratio, Q, from {3, 5, 7}, and the di¬ 
mensionless spin parameter, y, from {0.75,0.5,0}, where 
the spins are always prograde, i.e., parallel to the orbital 
angular momentum. Currently, neither the typical mass 
nor the typical spin of stellar-mass black holes is known 
from observations. Thus, we perform simulations system¬ 
atically adopting various values of them along with equa¬ 
tions of state to predict possible outcomes of binary merg¬ 
ers. Here, Q = 3, 5, and 7 correspond to Mbh = 4.O5M0, 
6 . 75 M 0 , and 9 . 45 M 0 , respectively. The low-mass black 
hole with « 4M0 is consistent with an observation of a 
black hole-Be star binary [SS], which could evolve into a 
black hole-neutron star binary, whereas the existence of 
a mass gap around 3-5Mq is frequently debated [HI HT] . 
The middle-mass, « 7Mq, and massive, « IOM 0 , black 
holes are safely expected to exist from observations of 
x-ray binaries [861187] . The spin parameter is even less 
constrained than the mass is [88j . and we simply take 
various values within our computational capabilities (see 
Ref. [31] for simulations of a near-extremal black hole- 
neutron star binary). We pay, however, less attention to 
high-mass and low-spin black holes. This is because such 
black holes are not able to disrupt companion neutron 
stars before they reach the innermost stable circular or¬ 
bit [H], and thus the merger process is essentially the 
same as that of binary black holes [53] . We also do not 
pay attention to retrograde spins, i.e., antiparallel to the 
orbital angular momentum, irrespective of the value of Q 
due to the same reason. 

Table |ll| also shows the adaptive-mesh-refinement grid 
structure for each simulation. We always choose {lc,lf) = 
(5,4) for Q = 3 and (4, 5) for Q = 5 and 7. In all the 
cases, the hydrodynamic evolution equations are solved 
only within L/2 « 1500 km for one side. Because it turns 
out later that a typical velocity of dynamical ejecta is 
0.2-0.3c, the ejecta motion can be safely tracked over 
~ 10 ms. At the same time, the box size is larger than 
the initial gravitational wavelength, and thus outgoing- 
wave boundary conditions are appropriate there as far 
as the gravitational wavelength is covered by > 10 grid 
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TABLE II. Key parameters of initial data and grid structures of simulations for models adopted in this study. The names of 
models represent the equation of state (EOS), mass ratio {Q), and dimensionless spin parameter of the black hole (x)- moQo, 
Mo, and Jo are the dimensionless initial orbital angular velocity, Arnowitt-Deser-Misner mass, and orbital angular momentum 
of the system, respectively. As for grid configurations, Ic and If are the numbers of coarser nonmoving grids and of a half of 
finer moving grids, respectively. The grid spacing at the finest level for A = 60 (fiducial resolution) is shown in physical units 
as well as a value normalized by mo « 2{Q + 1) km. The grid number assigned to the semimajor diameter of the neutron star 
is given by 7?diam/Aa; for the direction along the binary separation. The box size L is shown in physical units as well as a value 
normalized by the initial gravitational wavelength Ao = flo/rr. 


Model 

EOS 

Q 

X 

moflo 

Mo[Mq] 

MM@] 

Ic 

If 

Aa: (m) 

Axjmo 

-Rdiam/^^ 

L (km) 

L/Xo 

APR4-Q3a75 

APR4 

3 

0.75 

0.036 

5.35 

18.74 

5 

4 

162 

0.0203 

102 

2486 

3.6 

ALF2-Q3a75 

ALF2 

3 

0.75 

0.036 

5.35 

18.74 

5 

4 

186 

0.0233 

102 

2858 

4.1 

H4-Q3a75 

H4 

3 

0.75 

0.036 

5.35 

18.74 

5 

4 

209 

0.0263 

102 

3215 

4.6 

MSl-Q3a75 

MSI 

3 

0.75 

0.036 

5.35 

18.74 

5 

4 

228 

0.0286 

101 

3501 

5.0 

APR4-Q3a5 

APR4 

3 

0.5 

0.036 

5.35 

19.15 

5 

4 

162 

0.0203 

102 

2486 

3.6 

ALF2-Q3a5 

ALF2 

3 

0.5 

0.036 

5.35 

19.15 

5 

4 

186 

0.0233 

102 

2858 

4.1 

H4-Q3a5 

H4 

3 

0.5 

0.036 

5.35 

19.15 

5 

4 

209 

0.0263 

102 

3215 

4.6 

MSl-Q3a5 

MSI 

3 

0.5 

0.036 

5.35 

19.15 

5 

4 

228 

0.0286 

101 

3501 

5.0 

APR4-Q3aO 

APR4 

3 

0 

0.036 

5.35 

19.98 

5 

4 

162 

0.0203 

102 

2486 

3.6 

ALF2-Q3aO 

ALF2 

3 

0 

0.036 

5.35 

19.98 

5 

4 

186 

0.0233 

102 

2858 

4.1 

H4-Q3a0 

H4 

3 

0 

0.036 

5.35 

19.98 

5 

4 

209 

0.0263 

102 

3215 

4.6 

MSl-Q3aO 

MSI 

3 

0 

0.036 

5.35 

19.98 

5 

4 

228 

0.0286 

101 

3501 

5.0 

APR4-Q5a75 

APR4 

5 

0.75 

0.040 

8.04 

30.13 

4 

5 

158 

0.0132 

102 

2429 

2.6 

ALF2-Q5a75 

ALF2 

5 

0.75 

0.040 

8.04 

30.13 

4 

5 

181 

0.0152 

102 

2786 

3.0 

H4-Q5a75 

H4 

5 

0.75 

0.040 

8.04 

30.13 

4 

5 

205 

0.0171 

101 

3144 

3.3 

MSl-Q5a75 

MSI 

5 

0.75 

0.040 

8.04 

30.13 

4 

5 

219 

0.0183 

102 

3358 

3.6 

APR4-Q5a5 

APR4 

5 

0.5 

0.040 

8.04 

30.99 

4 

5 

158 

0.0132 

102 

2429 

2.6 

ALF2-Q5a5 

ALF2 

5 

0.5 

0.040 

8.04 

30.99 

4 

5 

181 

0.0152 

102 

2786 

3.0 

H4-Q5a5 

H4 

5 

0.5 

0.040 

8.04 

30.99 

4 

5 

205 

0.0171 

101 

3144 

3.3 

MSl-Q5a5 

MSI 

5 

0.5 

0.040 

8.04 

30.99 

4 

5 

219 

0.0183 

102 

3358 

3.6 

APR4-Q7a75 

APR4 

7 

0.75 

0.044 

10.73 

40.96 

4 

5 

153 

0.0096 

103 

2358 

2.1 

ALF2-Q7a75 

ALF2 

7 

0.75 

0.044 

10.73 

40.96 

4 

5 

179 

0.0112 

101 

2751 

2.4 

H4-Q7a75 

H4 

7 

0.75 

0.044 

10.73 

40.96 

4 

5 

200 

0.0125 

102 

3072 

2.7 

MSl-Q7a75 

MSI 

7 

0.75 

0.044 

10.73 

40.96 

4 

5 

215 

0.0135 

102 

3301 

2.9 

APR4-Q7a5 

APR4 

7 

0.5 

0.044 

10.73 

42.35 

4 

5 

154 

0.0097 

102 

2372 

2.1 

ALF2-Q7a5 

ALF2 

7 

0.5 

0.044 

10.73 

42.35 

4 

5 

179 

0.0112 

102 

2743 

2.4 

H4-Q7a5 

H4 

7 

0.5 

0.044 

10.74 

42.35 

4 

5 

201 

0.0126 

101 

3086 

2.7 

MSl-Q7a5 

MSI 

7 

0.5 

0.044 

10.74 

42.35 

4 

5 

217 

0.0136 

101 

3329 

2.9 


points. 


E. Diagnostics 

1. Ejecta 

We analyze global ejecta properties by integrals over 
unbound material [33]. We define the ejecta to be un¬ 
bound material identified by a criterion ut < —1, which 
becomes correct for a particle moving along its geodesics 
in a stationary spacetime. Because we are handling a 
fluid in a dynamical spacetime, this criterion is only ap¬ 
proximate and becomes especially poor in the vicinity 
of remnant black hole-disk systems. Our computational 
domains always extend to > 1000 km, where the grav¬ 
itational potential in geometrical units is < 0.01-0.02, 
and thus we expect that typical errors associated with 
this approximate criterion are a few percent. Strictly 
speaking, hut rather than Ut is a conserved quantity as¬ 
sociated with a fluid in a stationary spacetime. We check 


that the results depend only weakly on the choice of cri¬ 
teria, because shock heating does not play an important 
role in dynamical mass ejection from black hole-neutron 
star binaries (see Sec. IIIB). In consideration of the fact 
that our current simulations do not incorporate any pro¬ 
cess other than shocks responsible for heating and cooling 
such as neutrino interaction, we decide to neglect thermal 
effects for the purpose of classification. Because h > 1 by 
definition, our estimates should be regarded as conserva¬ 
tive. In addition, this allows us to compare our results 
directly with those of existing studies in numerical rela¬ 
tivity (e.g.. Refs. [3^l44] ). 


The rest mass outside the apparent horizon including 
both bound and unbound portions is computed by the 
integral 


Mr>rAH = / P*d^x, (14) 

J r>rAH 

where tah is the angle-dependent coordinate radius of 
the apparent horizon. The ejecta mass is defined by an 
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unbound portion of the rest mass as 


Mei = 


p* (fx. 


'r>T-AH,Ut<-l 


We also define the bound mass by 


(15) 


We also compute the linear momentum of ejecta, which 
indicates the degree of ejecta anisotropy. Components of 
the linear momentum of the ejecta may be defined by 

7^ej,2 — / p^Uid (^1) 

J r>rAH.'Ut< —1 


Mbd = / ptd^x 

• 7 r>rAH.“t>-l 

= iWr>rAH ~ -^ejj ( 16 ) 

which may be composed of the remnant disk and fall¬ 
back material. We do not, however, rigorously distin¬ 
guish these two components due to the absence of rea¬ 
sonable criteria. 

The kinetic energy of ejecta Tej is defined following 
Ref. [44]. First, the total energy of the ejecta is defined 
by 


where the z component vanishes identically due to the 
equatorial symmetry in this study. The magnitude of 
the linear momentum is given by 


— 



( 22 ) 


and the center-of-mass velocity of the ejecta may be de¬ 
fined by 


(23) 


F'ej = / p^ed^x, (17) 

J r>rAH.«t< —1 

whereas the gravitational binding energy is not (and can¬ 
not be in general relativity) appropriately subtracted. 
Next, the internal energy of the ejecta is defined by 

C/ej = / p*ed^x. (18) 

d r>rAii,u-t<-l 

Finally, the kinetic energy of the ejecta may be defined 
by subtracting the rest mass and internal energy from 
the total energy as 


which we call the bulk velocity in this paper. When 
the system is symmetric with respect to the equatorial 
plane, the bulk velocity vanishes if (but not only if) 
the ejecta is axisymmetric. A relation Uej < Uave al¬ 
ways holds. If the ejecta is modeled by an axisymmet¬ 
ric outflow truncated at an opening angle ipej, we have 
= Vave sin(vJej/2)/(<;3ej/2). Thesc measures suffer from 
the gravitational binding energy in the vicinity of black 
hole-disk systems as Tej and Uave do. Thus, they should 
also be estimated at a distant region. The propagation 
direction of the ejecta with respect to our coordinate sys¬ 
tem may be characterized by an angle defined from the 
linear momentum, 


Tej _ Tej Mej Tej- (19) 

Although the internal energy is likely to be converted to 
the kinetic energy in the long run, we do not count Tej as 
a part of Tej in this study. This does not affect the results, 
because Uej is smaller than Tej by orders of magnitude. 
Using the mass and kinetic energy of the ejecta, we may 
also define their average velocity as 


'^ave — 



( 20 ) 


using the Newtonian relation. It should be cautioned 
that the kinetic energy and average velocity defined in 
this manner are not calculated taking the gravitational 
binding energy associated with remnant black hole-disk 
systems into account. This implies that these measures 
overestimate asymptotic values when evaluated in the 
vicinity of black hole-disk systems independently of the 
validity of Ut < —1 and that they are reliable only for 
distant regions. For this reason, we typically measure 
the quantities of the ejecta at 10 ms after the onset of 
merger, when the dominant portion of the ejecta leaves 
the central region but still resides in our computational 
domains. 


d>ej = arctan [ j . (24) 

V ^ej,x J 

In addition to these integral quantities, the mass spec¬ 
trum with respect to the asymptotic velocity, or simply 
the velocity distribution of the ejecta, is estimated. The 
asymptotic velocity of each fluid element is defined from 
an asymptotic Lorentz factor —ut as 


V = Jl- 


{-utY ■ 


(25) 


Here, we use —Ut instead of the Lorentz factor seen from 
the Eulerian observer, au*, because the latter predicts 
the lower end of ejecta velocity to be the local escape 
velocity rather than zero. To derive the velocity distri¬ 
bution, we only analyze unbound material on the equa¬ 
torial plane and rescale the total mass to Mgj measured 


over the full region by Eq. (15). To compensate this ge¬ 


ometrical restriction, the mass of each fluid element is 
weighted by the distance from the coordinate origin, r, 
when computing the total mass of unbound material on 
the equatorial plane. This procedure is acceptable for 
black hole-neutron star binary mergers, because material 
ejected dynamically from neutron stars is concentrated 
around the equatorial plane (see Sec. IIIB). 









2. Fallback material 


We estimate fallback rates of bound material based 
on Newtonian relations [90]. The motion of the bound 
material is assumed to follow a ballistic trajectory de¬ 
termined by the energy and angular momentum of each 
fluid element. For this purpose, we only analyze bound 
material on the equatorial plane and rescale the total 
mass to Mbd measured over the full region by Eq. (16) 


in a similar manner to the computation of the velocity 
distribution of the ejecta. 

A fluid element on each grid point of the second-largest 
{I = 1) domain, which is the largest domain where the hy¬ 
drodynamic evolution equations are solved, is identified 
as an isolated test particle with the mass p{/S.xi)^ ne¬ 
glecting the spacetime curvature at a selected time slice. 
The specific energy (excluding the rest mass) E and spe¬ 
cific angular momentum J of the particle are estimated 
to be 


E = -ut - 1 , J = 


(26) 


where we only consider bound material identified by 
Ut > —1 and therefore E < 0. We neglect h — 1 in the 
same manner as the classification of the bound and un¬ 
bound material. The azimuthal velocity, u^p, is defined 
from Cartesian components by transformation with re¬ 
spect to the coordinate origin, which does not correspond 
exactly to the black-hole position nor center of mass (see 
the discussions in Sec. HIE). Assuming the presence of a 


central mass Me, the semimajor axis and eccentricity of 
the orbit are given by 


afb = - 


2E 


efb 



(27) 


in Newtonian gravity. Accordingly, the periapsis and 
apoapsis distances are given by Vp = afb(l — Cfb) and 
Ca = afb(l + efb), respectively. 

We define the fallback time of each particle to be the 
duration to reach the periapsis. The particle is assumed 
to obey the Newtonian equation of motion. 


dr Ur /„ ~ 2Mc 

= I—I \/ 2^^ H - 2 

at |Ur-| V ^ ’’ 


(28) 


regarding Ur as the radial velocity. This equation can 
be integrated analytically to give the fallback time for a 
particle at r = as 


Specifically, /(r^) is Pfb/4. It would be useful to recall 
that the orbital period is given by 


Pfb 




(31) 


Eor a particle with Cfb = 0, which appears in the central 
region, we simply set tfh = Tfb/2. Physically, compo¬ 
nents with efb ~ 0 should be regarded as the accretion 
disk rather than fallback material, while we do not have 
quantitative criteria to distinguish them. Such a par¬ 
ticle does not contribute in any way to the long-term 
fallback rate due to its short orbital period. We apply 
the same remedy for a particle that happens to satisfy 
< 0 and/or < Vp due to numerical errors, approxi¬ 
mate identification of the azimuthal velocity, or abuse of 
Newtonian relations. In this study, M^ is always approxi¬ 
mated by Too ignoring the energy loss due to gravitational 
waves and existence of the mass outside the black hole, 
Mr>rAH ■ We checked that the results depend only weakly 
on the precise value of Me- 

Einally, the fallback rate is computed by dividing the 
material into small segments according to the fallback 
time as 


Mfb(t) = 


AM{t) 
At{t) ’ 


(32) 


where AM{t) is the mass of the fluid elements satisfy¬ 
ing t < tfh < t + At and At{t) is arbitrarily chosen to 
be « t/IO. When we evaluate AM{t), the mass of each 
fluid element is weighted by r in the same way as done in 
the computation of the velocity distribution of the ejecta. 
It should be cautioned that, however, Mfb does not nec¬ 
essarily correspond to the black-hole accretion rate nor 
electromagnetic luminosity, because a part of the fallback 
material may be blown off from the disk as a wind or en¬ 
velope m- We do not discuss the fate of the fallback 
material in this study. 


3. Black hole 

Properties of remnant black holes are estimated by 
integrals on apparent horizons as in our previous work 
dSl [5nH55] . Assuming that the spacetime is approxi¬ 
mately stationary, the black-hole mass is estimated by 


ifb = ^ + ^[/(ra)-/(r.)], (29) 

2 \Ur\ 


MBH,f 


a 

47r ’ 


(33) 


where Pfb = is the orbital period and 


I(r) 


\/2Er'^ F 2Mrr - P 
~ 2E 


where Cg is the equatorial circumferential radius of the 
apparent horizon. The spin parameter of the remnant 
black hole Xi is estimated via the relation of Kerr black 
holes. 


Mg 

— aresm 

V-8E3 


/ 2Er + Mc\ 
y MgCfb j ' 


( 30 ) 


Cp /xf \ 

Ce tt \2f+)' 


(34) 


















9 


where Cp is the polar circumferential radius, f+ = l + 
x/i" — Xf is the normalized radius of the outer horizon, 
and E{z) is an elliptic integral defined by 



Comparisons among different estimates of the spin pa¬ 
rameter suggest that the systematic error associated with 
this method is Axf < 0.01 [29l [50l [51], and we do not 
repeat them here. 


4- Gravitational waves 


Our method to compute gravitational waves and re¬ 
lated quantities is summarized as follows (see Appendix 
B of Ref. [nH for the details). We extract the Weyl scalar 
'I '4 at « 4OOM0 from the coordinate origin by projecting 
onto spin-weighted spherical harmonics with i G {2,3,4} 
and extrapolate them to null infinity by a perturbative 
method |93j . The energy, linear momentum, and angu¬ 
lar momentum carried by gravitational waves are com¬ 
puted by integrating 4^4 in time |94j . The time integra¬ 
tion for calculating them and for deriving gravitational 
waveforms are performed by a fixed frequency integra¬ 
tion method [95] . Because we always impose the equa¬ 
torial symmetry, we only consider the 2 component for 
the radiated angular momentum and denote it as AJqw- 
The radiated linear momentum, which only has the x 
and y components, is decomposed into the magnitude 
APqw and angle 4 )gw in the same way as the ejecta [see 
Eqs. (221 and (24), respectively]. The radiated energy is 
denoted as AEgw- 


III. RESULT OF SIMULATIONS 

In this section, we present the results of our numerical 
simulations. Numerical values of characteristic quanti¬ 
ties are shown in Table m to which we refer repeatedly 
throughout this section. These values are estimated con¬ 
sistently at 10 ms after the onset of merger]^ 


A. Overview of merger dynamics 

We begin with a brief review of the dynamics of black 
hole-neutron star binary mergers (see Ref. [25| for de¬ 
tails) . Black hole-neutron star binaries evolve as a result 
of gravitational radiation reaction and eventually merge. 
Our initial conditions are chosen to evolve for ^ 3.5 to 


® We define the time of the onset of merger, tmerge, by the condi¬ 
tion that a part of neutron-star matter of mass O.OIMq falls into 
the apparent horizon in this and also previous work 129115211531 . 


7.5 orbits before the merger, where the exact numbers de¬ 
pend on model parameters. Eccentricities are estimated 
to be e ^ 0.01-0.02 for all the models using methods de¬ 
scribed in Ref. [S2], and they introduce uncertainties of 
the same order in the ejecta properties (see App. A 4 for 
the estimate). 

The fate of the system after the merger is determined 
primarily by competition between the orbital separation 
at which tidal disruption occurs (hereafter, the tidal dis¬ 
ruption radius), rtd, and the radius of the innermost sta¬ 
ble circular orbit, tisgo- If Usco is larger than rjd, no 
appreciable tidal disruption occurs, and the neutron star 
is simply swallowed by the black hole. In this case, the 
remnant disk, fallback material, and ejecta are all negli¬ 
gible for our astrophysical interest. Although we do not 
pay particular attention to such cases in this study, mod¬ 
els like APR4-Q3aO and ALF2-Q7a5 fall into this cate¬ 
gory (see the next paragraph). By contrast, if rtd is larger 
than risco j part of the disrupted material spreads around 
the black hole in the form of a tidal tail, while more than 
a half is still swallowed. Figure[^shows rest-mass density 
profiles on the equatorial plane in the central region at 
selected time slices for H4-Q5a75 as an example of this 
category. Material that remains outside the apparent 
horizon can be divided into bound and unbound mate¬ 
rial, and the former always dominates the latter for the 
models considered in this study[^ The bound material 
may be further divided into disk and fallback compo¬ 
nents. The unbound component is generated primarily 
by tidal torque exerted on the elongated neutron star dur¬ 
ing tidal disruption, and details of the dynamical mass 
ejection process are described separately in Sec. jlllB 


Appreciable tidal disruption occurs when (i) the 
neutron-star equation of state is stiff and the compact¬ 
ness is small, (ii) the mass ratio is small, and/or (hi) 
the black-hole spin is large (for a prograde orbit). These 
three conditions are reflected in our naming convention 
of the models. Note that, if we presume Mns to be fixed, 
condition i can be rephrased as “the neutron-star ra¬ 
dius is large” and condition ii as “the black-hole mass 
is small.” On one hand, rtd is expected to scale in the 
same way as the mass-shedding radius r^s, which is de¬ 
termined by the condition that the black-hole tidal force 
becomes equal to the neutron-star self-gravity at the stel¬ 
lar surface (see, e.g.. Ref. [2H]), 


rtd oc r^s 




NSj 


(36) 


and the dependence on the black-hole spin is not very 
strong [MlllTj- On the other hand, rigco is written as 
nsco(x)-^BH, where fisco(x) is a decreasing function 
of the dimensionless spin parameter, x [IH]- Recalling 
Rns/IWbh = 1 /(CQ), we expect the ratio to satisfy 


rtd _1_ 

nsco CQ’^/^nscoix)' 


(37) 


^ Hierarchy among the swallowed mass, bound mass, and unbound 
mass could change for extreme binary parameters m- 














10 


TABLE III. Characteristic physical quantities of the material measured at 10 ms after the onset of merger for our fiducial, 
N = 60 runs. Mr>rAii is the rest mass outside the apparent horizon. Mbd and Mej are the bound and unbound masses, 
respectively, and the unbound material is identihed as the ejecta. Note that Mr>rAH = Mbd + Mej. Tej and Pej are the 
kinetic energy and linear momentum of the ejecta, respectively. Uave = \/ 2Tej /Mej and Uej = Pej/Mej are the average and bulk 
velocities of the ejecta, respectively. 


Model 

^r>rjxH [-^o] 

Mbd[M0] 

Mej [Mq] 

Tej (erg) 

Tej[M0] 

't’ave 

Vej 

APR4-Q3a75 

0.19 

0.18 

0.01 

5 X 103*3 

2 X 10-3 

0.23 

0.19 

ALF2-Q3a75 

0.27 

0.23 

0.05 

3 X 10®3 

9 X 10-3 

0.25 

0.21 

H4-Q3a75 

0.33 

0.29 

0.05 

2 X 10®3 

9 X 10-3 

0.24 

0.20 

MSl-Q3a75 

0.35 

0.28 

0.07 

4 X 10®3 

0.01 

0.25 

0.21 

APR4-Q3a5 

0.08 

0.08 

2 X 10“^ 

1 X 103*3 

4 X 10-3 

0.21 

0.17 

ALF2-Q3a5 

0.19 

0.17 

0.02 

1 X 10*3 

5 X 10-3 

0.24 

0.20 

H4-Q3a5 

0.24 

0.21 

0.03 

1 X 10®3 

6 X 10-3 

0.23 

0.20 

MSl-Q3a5 

0.26 

0.21 

0.05 

3 X 10*3 

0.01 

0.24 

0.21 

APR4-Q3aO 

4 X 10”“ 

4 X 10"^ 

2 X 10"'’ 

6 X 10'‘''' 

1 X 10-3 

0.20 

0.08 

ALF2-Q3aO 

0.03 

0.03 

3 X 10-3 

1 X 10**3 

3 X 10-3 

0.22 

0.11 

H4-Q3a0 

0.10 

0.10 

6 X 10-3 

3 X 10**3 

1 X 10-3 

0.22 

0.18 

MSl-Q3aO 

0.16 

0.14 

0.02 

8 X 10**3 

3 X 10-3 

0.23 

0.19 

APR4-Q5a75 

0.07 

0.06 

8 X 10-3 

5 X 10**3 

8 X 10-3 

0.25 

0.10 

ALF2-Q5a75 

0.24 

0.20 

0.05 

3 X 10*3 

0.01 

0.28 

0.21 

H4-Q5a75 

0.32 

0.27 

0.05 

3 X 10*3 

0.01 

0.27 

0.22 

MSl-Q5a75 

0.36 

0.28 

0.08 

6 X 10*3 

0.02 

0.28 

0.23 

APR4-Q5a5 

5 X 10"^ 

4 X 10"^ 

9 X 10-3 

4 X 10'33 

5 X 10-3 

0.23 

0.05 

ALF2-Q5a5 

0.04 

0.03 

0.01 

8 X 10**3 

7 X 10-3 

0.27 

0.06 

H4-Q5a5 

0.14 

0.12 

0.02 

1 X 10*3 

4 X 10-3 

0.26 

0.19 

MSl-Q5a5 

0.23 

0.18 

0.05 

3 X 10*3 

0.01 

0.27 

0.21 

APR4-Q7a75 

2 X 10"=* 

2 X 10"^ 

5 X 10-^ 

4 X 10'‘'3 

3 X 10-3 

0.27 

0.06 

ALF2-Q7a75 

0.07 

0.05 

0.02 

2 X 10*3 

2 X 10-3 

0.29 

0.07 

H4-Q7a75 

0.19 

0.16 

0.04 

3 X 10*3 

7 X 10-3 

0.29 

0.19 

MSl-Q7a75 

0.30 

0.23 

0.07 

5 X 10*3 

1 X 10-3 

0.30 

0.23 

APR4-Q7a5 

1 X 10“'’ 

1 X 10"'’ 

3 X 10-'= 

1 X lO'** 

1 X 10-' 

0.23 

0.04 

ALF2-Q7a5 

5 X 10""^ 

3 X lO"'^ 

2 X 10-"^ 

1 X 10'‘® 

9 X 10-® 

0.27 

0.05 

H4-Q7a5 

6 X 10"^ 

3 X 10"® 

3 X 10-3 

3 X 10**3 

2 X 10-3 

0.29 

0.06 

MSl-Q7a5 

0.04 

0.02 

0.02 

1 X 10*3 

1 X 10-3 

0.30 

0.07 



FIG. 1. Rest-mass density profile on the equatorial plane in the central region at selected time slices for H4-Q5a75. Black filled 
circles show the interior of apparent horizons. The left, middle, and right panels correspond to a late inspiral phase, tidal tail 
formation, and quasistationary remnant accretion disk, respectively. The merger sets in at t = 26.49 ms for this model. 


and a large value of this ratio should signal apprecia¬ 
ble tidal disruption. This expectation has been verihed 
by previous studies of disk formation and gravitational- 
wave emission [28) . and Table III indicates that dynami¬ 
cal mass ejection also becomes efficient when these three 


parameters {C, Q, and x) are advantageous for tidal dis¬ 
ruption. The dependence of the ejecta properties on these 
parameters is discussed in more detail in Sec. [Ill C| 
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B. Mass ejection process and morphology 


We first explain mechanisms of dynamical mass ejec¬ 
tion and general properties of the ejecta by closely inves¬ 
tigating APR4-Q3a75 in Sec. |III BT| Mass ejection mech¬ 
anisms and qualitative trends are the same for all the 
black hole-neutron star binary models simulated in this 
study, whereas differences in (semi)quantitative proper¬ 
ties are found. We next discuss differences in the ejecta 
geometry among models in Sec. IIIB2 Characteristic 


quantities and their differences are described in Sec. Ill C 


1. Case study: APR4-Q3a75 

Figure [^depicts the typical process of dynamical mass 
ejection at tidal disruption for model APR4-Q3a75. Once 
tidal disruption sets in, the neutron star is drastically 
elongated and forms a tidal tail. While the high-density 
innermost part is immediately swallowed by the black 
hole, the outer part spreads to a distant orbit and lags 
behind. Thus, the tidal tail exhibits a trailing one-armed 
spiral structure, and the black hole exerts tidal torque on 
the tail, increasing its orbital angular momentum. The 
outer part of the tail moves further outward due to the 
gain of the angular momentum, and the outermost part 
acquires enough kinetic energy to become unbound from 
the system, as marked by black curves in Fig. In the 
course of this process, the pressure gradient in the tidal 
tail should also boost the outer part. This angular mo¬ 
mentum transport proceeds in an unstable manner until 
the tidal tail winds around the black hole and collides 
with itself to form a nearly axisymmetric black hole-disk 
system. This mechanism generates most of the dynami¬ 
cal ejecta as well as bound material which eventually falls 
back to the black hole-disk system. 

Although a small amount of unbound material appears 
to be ejected toward ip ^ —45° with a large velocity 
in Fig. where p is the azimuthal angle in spherical 
coordinates, this appears to be an artifact created by 
the artificial atmosphere and finite grid resolutions as we 
discuss in Appendix |A 5| This observation is consistent 
with Ref. |33j . The mass, energy, and linear momentum 
of this component is negligible compared with the main 
component discussed in the previous paragraph, and thus 
the values shown in Table m are not affected. 

Dynamical mass ejection from black hole-neutron star 
binaries is anisotropic [49]. Figure shows that the 
ejected material takes a crescentlike shape on the equato¬ 
rial plane during its early evolution for APR4-Q3a75. Al¬ 
though the relative size of the central region occupied by 
bound material will become negligible as the rear velocity 
approaches zero (see below), the ejecta never sweep out 
the whole equatorial plane. Furthermore, it is concen¬ 
trated around the equatorial plane and does not extend 
above the central black hole, because this mass ejection is 
driven by the tidal torque, which works most efficiently in 
the equatorial plane. This situation should be contrasted 


with dynamical mass ejection from binary neutron stars, 
in which quasiradial oscillations of remnant massive neu¬ 
tron stars eject an appreciable amount of material toward 
polar regions via shock interaction [Ji] . To elucidate the 
difference, we show the thermal part of specific internal 
energy, eth, in Fig. i As shocks do not play a role, 
the tidal tail including the ejecta is not heated signifi¬ 
cantly except for the self-colliding region of the tidal tail. 
The self-colliding shock interaction eventually thermal- 
izes and circularizes material in the central region, and a 
hot accretion disk is formed. We will discuss properties 
of the accretion disk later (see also Sec. IIIFl. Apparent 
heating at the outermost part of the tidal tail is caused 
by the artificial atmosphere and is thus spurious. 

Figures and [^suggest that the dynamical ejecta orig¬ 
inates from the outer core and crust of the neutron star 
retaining its very low electron fraction (the number of 
electrons per baryon). Ye < 0.1, at zero temperature 
[55] Because Mej for APR4-Q3a75 is comparable to 
the typical mass of neutron-star crusts, O.OIMq (see, 
e.g.. Ref. [99]), the ejecta stripped from the outermost 
part of the tidal tail in a highly nonspherical manner 
stems not only from the crust but also from the core 
(see also Fig. 3 of Ref. [31]). In fact, Mej for other bi¬ 
nary models can easily exceed the typical crust mass. 
Nevertheless, the ejecta would not come from the in¬ 
ner core, because the densest part of the neutron star 
is swallowed by the black hole and bound material sep¬ 
arates the black hole and ejecta. Thus, the dynamical 
ejecta should come mainly from the outer core and partly 
from the crust. The absence of shocks suggests that the 
low electron fraction of the outer core is not modified 
very much during dynamical mass ejection, and this is 
consistent with results obtained by previous smoothed- 
particle-hydrodynamics simulations [33H31]- Such ejecta 
are expected to be a promising site of r-process nucle¬ 
osynthesis producing predominantly second- and third- 
peak elements via fission cycling, while the production 
of first-peak elements may not accompany [1001 llOlj . It 
has to be cautioned that this estimation is speculative to 
some extent, because our simulations are performed with¬ 
out taking the electron fraction into account. We plan 
to revisit this topic with more sophisticated equations of 
state and neutrino transport schemes [102] . 

Figure]^ shows the long-term evolution of the dynam¬ 
ical ejecta in the distant region. This figure shows that 
the outer edge of the ejecta expands in a nearly homoge¬ 
neous manner after the angular momentum transport by 
the tidal torque ceases. The azimuthal component of ve¬ 
locity decreases approximately as r~^ due to angular mo¬ 
mentum conservation and soon becomes negligible com¬ 
pared to the radial component as shown in Fig. [^ This 


Identifying the origin of postmerger material is much more 
difficult in mesh-based simulations than in smoothed-particle- 
hydrodynamics simulations. Rigorous confirmation would re¬ 
quire postprocess calculations using Lagrangian tracer particles. 
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FIG. 2. Rest-mass density profile for APR4-Q3a75 during the tidal disruption and dynamical mass ejection. Note the different 
spatial and density scales among the columns. Black filled circles show the interior of the apparent horizons. Black arrows 
show the spatial component of the covariant four-velocity, m. Unbound components satisfying ut < —1 are marked by black 
curves, and we checked that contours marking hut < — 1 nearly overlap with them. The top, middle, and bottom rows are the 
xy, xz, and yz planes, respectively. The merger sets in at t = 18.48 ms for this model. This figure should be compared with 
Figs. 3-5 of Ref. |44) where binary neutron stars are studied. 


implies that the kinetic energy of the ejecta is dominated 
by the radial velocity, and thus the average velocity, Uave, 
estimated from the kinetic energy approximately equals 
the typical radial velocity. Opening angles of ejecta in 
the equatorial and also meridional (not shown in Fig. 
planes are approximately conserved, because the direc¬ 
tion of velocity does not change appreciably once hydro- 
dynamic interaction becomes negligible. Note that en¬ 
ergy injection by the r-process heating will moderately 
change the ejecta geometry m- 

Figure also shows that the radial thickness of the 
dynamical ejecta increases in the long-term evolution of 
the ejecta, because the ejecta head is faster than the 
rear. Specifically, the head will maintain a velocity on 
the order of the escape velocity of the neutron star, while 
the rear velocity will approach zero (separation of bound 
and unbound components) as the material climbs up the 
gravitational potential well. The radius of the central 
bound region will become negligible compared to the ra¬ 
dial thickness of the dynamical ejecta for exactly the 
same reason. 

After disk formation, unbound material is newly gen¬ 


erated from the disk region due to shock heating. Figure 
shows the shock heating-driven disk outflow on the xz, 
meridional plane. When the tidal tail collides with itself, 
shock interaction increases Sth near the contact surface. 
Because the rest-mass density is not high in the relevant 
region, thermal pressure dominates the cold-part pres¬ 
sure]^ The heated material expands, and some material 
is puffed up off the equatorial plane. In addition, shock 
interaction circularizes incoming tail material, and thus 
the disk region extends radially. Cold fallback material 
eventually accumulates and circulates around the outer 
edge of the hot disk material, as is visible from the right 
panel of Fig. at a; « 120 km. When the accumulated 
cold material becomes very massive, shocks develop be¬ 
tween the cold and hot material. Shock heating occurs 
continually at the outer edge of the disk due to this in¬ 
teraction, and material is also puffed up there. Material 


® This should correspond to the dominance of gas and radiation 
pressure over electron degeneracy pressure. 
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FIG. 3. Profile of the thermal part of specific internal en¬ 
ergy, £th, for APR4-Q3a75. The time slice is taken to be 
the same as the right column of Fig. The black filled cir¬ 
cle, arrows, and curves have the same meanings as those in 
Fig. i. Green dashed curves show contours of p = 10®, 10^°, 
and 10^^ gcm“®. This figure should be compared with Figs. 7 
and 8 of Ref. |44| where binary neutron stars are studied, tak¬ 
ing different spatial scales into account. 


off the equatorial plane exhibits (seemingly) random mo¬ 
tion, and a part of it collides with another part. Finally, 
some of the material is ejected from the system as hot 
blobs, and the rest eventually falls back to the disk sur¬ 
face. In contrast to dynamical mass ejection due to the 
tidal torque, this mechanism ejects material mainly to¬ 
ward nonequatorial (vertical) directions. As is evident 
from the left panel of Fig. however, the mass of the 
ejecta generated by this heating is much smaller than 


that by the tidal torque. The situation will change if 
magnetic fields [lOdj . neutrino heating [Ml 1105) . and/or 
nuclear interactions [1061 I107j are taken into account, 
whereas a significant fraction of the disk material has to 
be ejected to dominate over dynamical mass ejection. 


8. Variety of ejecta morphology 

The ejecta geometry may be characterized by an open¬ 
ing angle in the equatorial plane, (pej, and that in the 
meridional plane, dej, where the latter is defined to refer 
only to material with z > 0, taking the equatorial sym¬ 
metry into account. In the nearly spherical mass ejection 
expected for supernovae and binary neutron star mergers 
HU, ipej and 0ej should be regarded as 27r and 7r/2, respec¬ 
tively. We give estimates based on analytic arguments of 
the opening angles for black hole-neutron star binaries in 
Appendi x |B] to compare with numerical results. 

Figure [6[shows the morphology of the dynamical ejecta 
on the equatorial plane for various models. This figure 
implies that a softer equation of state, a larger mass ra¬ 
tio, and a smaller black-hole spin lead to a larger value 
of (pej when other parameters are fixedj^ In particular 
for the case in which mass ejection is not very substan¬ 
tial, an unbound portion revolves more than one orbit 
(tpej > 27r) taking a spiral shape at generation, and rear- 
end collisions occur in overlapping directions to form a 
ring shape. Traces of the rear-end collisions are observed 
as bumps on boundaries between bound and unbound 
material (black closed curves) in Fig. In these cases, 
the bulk velocity, Uej, is lower than O.Ic and is less than 
half of the average velocity, Uave (see Table |nl| ). The rea¬ 
son for this is that the ejecta linear momentum, Pej, is 
very small for nearly axisymmetric mass ejection. 

This catalog suggests that tpej tends to become large 
when tidal disruption occurs only weakly. This tendency 
does not agree with the estimate obtained by time-scale 
arguments in Appendix]^ A possible explanation of this 
tendency is the periastron advance in general relativity, 
which is pronounced when tidal disruption occurs very 
close to the innermost stable circular orbit m- As an 
extreme example, orbital parameters of a test particle can 
be finely tuned so that it experiences an arbitrarily large 
number of revolutions traveling near marginally stable 
orbits [TMfTTU] . Although the ejecta material cannot 
be finely tuned due to its finite spatial extent and does 
not experience infinitely many revolutions, i.e., (pej will 
not diverge, the dynamical ejecta should be able to have a 
large value of tpej if the mass ejection takes place near the 
innermost stable circular orbit. Indeed, tidal disruption 
should have occurred very close to the innermost stable 
circular orbit, i.e., rtd ~ nscOi when the ejecta mass 


APR4-Q3aO might seem to have a smaller than ALF2-Q3aO, 
but this simply reflects the fact that the ejecta of APR-Q3aO is 
extremely tiny. 
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FIG. 4. The same as Fig. but on the equatorial plane at late times in the distant region. The black curves and arrows have 
the same meanings as those in Fig.In this model, the ejecta linear momentum points toward <l>ej ~ —100°, i.e., close to the 
— y direction [see Eq. (|24[) for the definition]. 
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FIG. 5. Profile of rest-mass density (left) and thermal part of the specific internal energy (right) on the xz plane for APR4- 
Q3a75 after the disk formation. The black arrows and curves have the same meanings as those in Fig.|^ Green dashed curves 
show contours of p = 10®, 10^°, and 10^^ gcm“®. 


is small but nonnegligible. This is consistent with the 
tendency observed in Fig. 

From the observational viewpoint, dynamical ejecta 
with a large opening angle, tpej ^ 27r, may not be very im¬ 
portant, because a large opening angle is attained by the 
ejecta with a small mass, for which electromagnetic ra¬ 
diation is expected to be weak. Strong electromagnetic 
radiation should accompany substantial mass ejection, 
say Mej > O.OIMq, where fpej takes a value close to tt in 
most cases. However, substantial but nearly axisymmet- 
ric dynamical mass ejection such as that for ALF2-Q7a75 
is not completely excluded. 

The opening angle in the meridional plane does not dif¬ 
fer very much among models as far as substantial mass 
ejection occurs. Figure[7]shows the morphology of the dy¬ 
namical ejecta on the meridional plane for various mod¬ 
els. This figure shows that the opening angle 0ej takes 
values between 1/5 and 1/3 (or 10° and 20°) for cases 
with Mej > O.OIMq. The variation of dej up to a fac¬ 
tor of < 2 is observed among models with substantial 
mass ejection, but the ejecta driven by the tidal torque 
never extend to, say, 0ej > 30°. At the same time, 0ej 


is very small when mass ejection is not efficient. Hence, 
sphericity is never achieved even approximately for cases 
considered in this study. This figure also suggests that 0ej 
tends to become small when Q is large. This is consistent 
with the analytic expectation presented in Appendix [B| 


C. Characteristic quantities of ejecta 


Here we discuss characteristic quantities of dynamical 
ejecta such as the mass and velocity, focusing on their 
dependence on binary parameters. As described in the 
beginning of this section, we measure ejecta quantities at 
10 ms after the onset of merger. To check that estima¬ 
tion at that time gives acceptable results, we first investi¬ 
gate time evolution of the ejecta quantities in Sec. ImcT] 


Next, we discuss the dependence in Sec. HIC2 
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FIG. 6. Rest-mass density profile on the equatorial plane for various models at « 10 ms after the onset of merger. Unbound 
components satisfying ut < — 1 are marked by black curves. From left to right, the equation of state is APR4, ALF2, H4, 
and MSI. From top to bottom, {Q,x) is (3,0), (3,0.5), (3,0.75), (5,0.75), and (7,0.75). The left panel on the third row is 
APR4-Q3a75 described in detail in Sec. |IIIB 1| Traces of rear-end collisions are found as bumps on inner black closed curves 
But tp ~ -K 12 for ALF2-Q3a0 and APR4-Q5a75 and at ip « 0 for APR4Q7a75 and ALF2Q7a75. 


1. Time evolution 

Figure!^ shows the time evolution of M^j, Tgj, and Pej 
for selects models. All these values suddenly increase 
right after the onset of merger. The time evolution in¬ 
dicates that most of dynamical mass ejection progresses 
over « 2 ms and that the evolution relaxes afterward ir¬ 
respective of the models. 


The ejecta mass settles to a quasistationary value 
within ~ 5 ms. This confirms the observation in 

Sec. |IIIB 1| that mass ejection due to disk activity does 
not contribute significantly to the total mass of the ejecta 
in our simulations. Therefore, the measurement of Mgj 
at 10 ms after the onset of merger is safely justified. 

The kinetic energy and linear momentum peak at 1- 
2 ms after the onset of merger and decrease afterward. 



































































16 


1200 
i 600 
0 

-1200 -600 0 600 1200 
y (km) 

1200 


—1—1—1—1—I—1—1—1—1—1—1— 

0.5c- APR4-Q3aO 

■ 12 1200 
® 11 

- 10 - 
! 8 I 600 

0.5c- ALF2-Q3a0 

■ 12 1200 
- 11 

- 10 .p 

! 8 *600 

0.5c - H4-Q3a0 

■ 12 1200 

* 11 

- 10 -p 

! 8 1. 600 

0.5c- MSl-Q3a0 ! 


-7 

■ 6 

J 5 0 


-7 

■ 6 

5 0 

. . 1 r-^ilT^lliCTl 

-7 ^ 

■ 6 

^5 0 

■ 1 ir 


-1200 -600 0 600 1200 
y (km) 


-1200 -600 0 600 1200 
y (km) 


-1200 -600 0 600 1200 
X (km) 


I 600 


0 


0 . 5 c - APR 4 - Q 3 a 5 




1200 


0 . 5 c - ALF 2 - Q 3 a 5 

600 - 



■ 


1200 r 


- 11 




- m 



- 

1 9 

- 8 

600 - 


- 7 

N 



■ 6 

5 


0 - 


0 . 5 c - H 4 - Q 3 a 5 


1200 


I 600 


MSl-Q3a5 


' ' 1 ' ‘ 1 ' ' 1 ' ' 

• 0.5c- APR4-Q3a75 ' 

■ 12 1200 
■ 11 

- 10 ^ 

—1—1—1—1—^—|—1—1—1—1—1— 

0.5c- ALF2-Q3a75 

■ 12 1200 
- 11 
- 10 

0.5c- H4-Q3a75 ' 

■ 12 1200 
■ 11 

- 10 -c? 

0.5c- MSl-Q3a75 ' ! 

- 

! 1 1 600 

- 

! 8 * 600 

- 

*8 600 

- ! 


-7 ^ 

■ 6 

J 5 0 


-7 

*5 0 


-7 « 

*5 0 



-1200 -600 0 600 1200 

y (km) 

1200 
I 600 

0 

-I 

1200 
I 600 
0 

-1200 -600 0 600 1200 

y (km) 


-1200 -600 0 600 1200 
y (km) 


-1200 -600 0 600 1200 
y (km) 


-1200 -600 0 600 1200 
X (km) 


-1200 -600 0 600 1200 
y (km) 


-1200 -600 0 600 1200 
y(km) 


-1200 -600 0 600 1200 
y (km) 


-1200 -600 0 600 1200 
X (km) 



1200 


600 - 
0 


-T—I—I—I— I -1—I—I—I—I—r- 

0 . 5 c - ALF 2 - Q 5 a 75 





■ 1^ 


1200 r 


S 11 




- 1(1 

I 


- 

1 9 

- 8 

600 - 


- 7 

EM 



■ 6 

5 


0 - 


0 . 5 c - H 4 - Q 5 a 75 


1200 


I 600 


MSl-Q5a75 




-1200 -600 0 600 1200 
X (km) 


-1200 -600 0 600 1200 
X (km) 


-1200 -600 0 600 1200 
y(km) 


—1—1—1—^—1—1—1—1—1—^—1— 
0.5c- APR4-Q7a75 ' 

■ 12 1200 
■ 11 

- 10 ^ 

! 1 1 600 

0.5c- ALF2-Q7a75 

■ 12 1200 
- 11 

- 10 - 
! ^ 1 600 

0.5c- H4-Q7a75 ' 

m 12 1200 

■ 11 

\ 1 g 600 

0.5c- MSl-Q7a75 I 


-7 

■ 6 

^5 0 


-7 ^ 

"5 0 


-7 ^ 

"5 0 



-1200 -600 0 600 1200 -1200 -600 0 600 1200 -1200 -600 0 600 1200 -1200 -600 0 600 1200 
X (km) X (km) y (km) x (km) 


FIG. 7. The same as Fig. but on a meridional, xz or yz plane, chosen to be the one which is closer to the direction of $ej. 
Unbound components satisfying ut < —1 are marked by black curves. From left to right, the equation of state is APR4, ALF2, 
H4, and MSI. From top to bottom, {Q,x) is (3,0), (3,0.5), (3,0.75), (5,0.75), and (7,0.75). The left panel on the third row is 
APR4-Q3a75 described in detail in Sec. |III B 1[ 


The reason of this decrease is that the ejecta lose en¬ 
ergy in climbing up the gravitational potential well of the 
central black hole-disk system. The Newtonian formulas 
indicate that Tej measured at 10 ms after the onset of 
merger overestimates its final value by (moMej/r)/Tej ^ 
30%-40% for models shown in Fig. and this is consis¬ 
tent with the later evolution. This will result in ^ 15%- 
20 % overestimation of the ejecta velocity, and thus this 
error has to be kept in mind in the following discussions, 
along with those described in Appendix]^ If we measure 
these values at < 5 ms after the onset of merger and use 
them as proxies for their final values, final ejecta veloci¬ 
ties can be overestimated nearly by 100%. Hence, a large 
computational domain is a prerequisite for an accurate 
study of mass ejection^] 


The amount of error depends on estimation methods. For exam¬ 
ple, the kinetic energy can also be defined by f pt{—ut — l)tfix 
(F. Foucart, private communication). 


2. Dependence on binary parameters 


We start by looking at the total mass remaining out¬ 
side the apparent horizon, Mr>rAH = -^bd+fHej, to check 
consistency with previous work. Figure plots Mry^AH 
measured at 10 ms after the onset of merger (presented in 
Table III) as a function of the compactness, C. This fig¬ 
ure supports the discussion in Sec. Ill A That is, a small 
neutron-star compactness, small mass ratio, and large 
black-hole spin increase the strength of the tidal disrup¬ 
tion resulting in the increase of Mr-yrAu- present 

simulations reproduce quantitatively the results of our 
previous simulations [29l [52l [53], as well as those by 
other authors (see Ref. [Ill] for a compilation). The de¬ 
pendence of Mr>rAH ^ approximately linear within 
the range studied here, until it levels off at < O.OIMq. 
This suggests that the effect of neutron-star properties 
on Mj-yrAH is reasonably captured by the compactness, 
C. 


The ejecta mass, Mgj, is correlated with the strength 
of the tidal disruption as i®! iii® depen¬ 

dence of Mej on binary parameters is complicated. Fig- 
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FIG. 8. Time evolution of the mass (top), kinetic energy 
(middle), and linear momentum (bottom) of the ejecta for 
selected models. 


To] shows Mej as a function of C. Plots of Tej and 


ure 

Pej exhibit similar behavior. For fixed values of Q and 
X, Mej increases as C decreases. This is qualitatively 
the same as supports the expectation that 

strong tidal disruption is accompanied by efficient mass 
ejection. However, the correlation is weaker between Mej 
and C than between Mr>rAH ^ for fixed values of 
Q and x- This suggests that the boundary separating 
bound and unbound material, Ut = —1, is not deter- 


FIG. 9. Mass remaining outside the apparent horizon com¬ 
posed of bound and unbound material, Mr>r ah ~ Mbd + Mej, 
as functions of the compactness, C. 
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FIG. 10. Ejecta mass, Mgj, as a function of the compactness, 

C. 


mined solely by the compactness but is also sensitive to 
the stellar structure. This observation is consistent with 
Ref. [33j . which found a similar fact by comparing their 
results with some of our results reported in Ref. m- 
It is reasonable that detailed properties of the equation 
of state could play an important role during dynamical 
mass ejection via effects such as the pressure gradient 
and/or central condensation. 

The ejecta mass, Mej, does not depend monotonically 
on the mass ratio, Q, for fixed values of C and x (see 
Fig. 10). The reason for this is that the ejecta tends to 
comprise a large fraction of material remaining outside 
the apparent horizon for a large mass ratio, especially 
when tidal disruption is weak and Mr>rAH is not very 
large. Figure [TT] shows the correlation between the ejecta 
mass, Mej, and bound mass, Mbd. This figure indicates 
that Mej does not decrease very rapidly with the decrease 
of Mbd (and equivalently Mr>r ah) ^ large value of 
Q. Specifically, Mej > O.OIMq can be achieved when 
Mbd ^ O.OIMq for Q = 7, while it is possible only when 
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FIG. 11. Correlation between the ejecta mass, Mej, and 
bonnd mass, Mbd. We restrict the region to Mej > O.OOIM© 
and Mbd > O.OIM©, where the results are astrophysically 
interesting. 


hole, where a half of the star is expected to become un¬ 
bound |113] . Although the qualitative trend is the same, 
dynamical processes should play a crucial role in realizing 
this dependence in the inspiral of black hole-neutron star 
binaries, because all the neutron-star material is bound 
to the system at the onset of tidal disruption. Note that 
the systematic error in Uave associated with the resid¬ 
ual gravitational binding described in Sec. EH C 1| is not 
likely to modify this tendency qualitatively, because all 
the values of v^ve are systematically overestimated. 

The dependence of the ejecta mass, M^j, on the black- 
hole spin, X, is simpler than that on C and Q (see Fig. 101. 
Namely, a large black-hole spin increases the amount of 
ejecta for fixed values of C and Q. We do not find signifi¬ 
cant dependence of Mej/Mbd on %. The average velocity, 
Vave, tends to increase as x increases. 

The ejecta mass, Mej, is correlated with the mass re¬ 
maining outside the apparent horizon, Mr^rAm indi¬ 
cated in Fig. El Quantitatively, we obtain 
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FIG. 12. Average velocity of the ejecta, Wave, defined from the 
kinetic energy, Tej, as a function of the compactness, C. 


Mhd ^ O.IMq for (5 = 3. The fact that mass ejection can 
be substantial even if tidal disruption is not very strong 
for a large value of Q is encouraging for electromagnetic 
counterpart searches, because astrophysical black holes 
are expected to prefer large mass ratios [86l |87j . 

The increase of Mgj/Mbd with the mass ratio, Q, im¬ 
plies that material remaining outside the apparent hori¬ 
zon tends to become more energetic when Q is larger. 
This speculation is supported by the fact that the av¬ 
erage velocity of the ejecta, r^ave, is larger for a larger 
mass ratio. Figure [T^ shows Uave as a function of C. A 
typical value of Uave is 0.22-0.25c for (5 = 3, and this 
rises to 0.25-0.28c for (5 = 5 and 0.28-0.3c for Q = 7. 
This can be ascribed to the higher energy of material re¬ 
maining outside the apparent horizon for a larger value 
of Q. The effect of the mass ratio on ejecta velocities 
via a gravitational potential is pointed out in the context 
of the tidal disruption of a main sequence star during 
a nearly parabolic encounter with a supermassive black 


Mej 

Mq 


(0.27 ±0.07) 


f ^r>rAH 

V Mq 


1.3±0.2 


(38) 


by fitting all the data shown in Table |III| with equal 
weights, where the range indicates the I-cr asymptotic 
standard error. If we fit the data of models with differ¬ 
ent values of Q separately, relations become 


Mej 

Me 


(0.41 ±0.14) {Mr^rAu/Me) 


1.8±0.3 


^ = <((0.23 ±0.06) {Mr>rAjMe) 
(0.15 ±0.02) iMr>rAH/Me) 


1.1±0.2 

0.73±0.09 


(Q = 3) 
(Q = 5) . 
(Q = 7) 
(39) 

It is evident that the power-law index is smaller for a 
larger value of Q, and thus the separate fitting may be 
more appropriate. These relations give us an approxi¬ 
mate estimate of Mej combined with a fitting formula for 
Mr>rAH provided in Ref. |111) . Sources of the error come 
from both simulations and fitting procedures, and only 
the latter is taken into account in Eqs. (38) and (39). 


D. Ejecta and envelope structure 

First in Sec. IIIIDll we investigate matter prohles on 
the equatorial plane, where most of the material resides. 
It includes disk, fallback, and ejecta components. Next, 
material distribution along the 0 axis is investigated in 
Sec. |IIID 2[ It will be important for gamma-ray bursts, 
because a hypothetical jet (or fireball) can achieve an 
ultrarelativistic velocity only if the baryon load is not 
very high [mi- Finally, we investigate the velocity dis¬ 
tribution of dynamical ejecta in Sec. IIIID31 which is 
required to predict electromagnetic radiation quantita¬ 
tively [71 ng. Detailed structures of material obtained 
from our simulations are not expected to be very real¬ 
istic, because the equation of state in a relevant regime 
is composed of a single zero-temperature polytrope and 
ideal-gas-like thermal correction. We still believe that 
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our results capture qualitative properties of the mate¬ 
rial structure, particularly for ejecta in distant regions 
where hydrodynamic interaction does not play an impor¬ 
tant role. 


1. Equatorial plane 

Figure shows density profiles along the x and y axes 
at 10 ms after the onset of merger for selected models. 
Corresponding snapshots are given in Fig. The mate¬ 
rial at r < 100 km is in an approximately axisymmetric 
state for all the models. This implies that accretion disks 
are formed in the central regions at this time. For a given 
value of rest-mass density in the disk region 

is higher when Q is smaller. The reason for this is that 
characteristic length scales are proportional to the total 
mass of the system, and thus to Q -I- 1. Accordingly, 
characteristic rest-mass density should be proportional 
to (Q + 1)”^ for a given value of Mr^rAH- This tendency 
was already reported in Ref. |29j. 

Density profiles outside the disk region depend signif¬ 
icantly on the azimuthal angle. On one hand, the rest- 
mass density steeply decreases along directions with no 
ejecta. In Fig.[^ the +y direction of APR4-Q3a75 and 
—X direction of H4-Q5a75 fall into this category. The +x 
direction of MSl-Q7a75 also corresponds to this case, but 
a high-density region is still observed up to « 500 km, be¬ 
cause the tidal tail has not fallen back and collided with 
itself yet in this direction. On the other hand, approxi¬ 
mately constant density plateaus extend up to ^ 1000 km 
along directions that the ejecta sweep. For example, the 
—X and —y directions of APR4-Q3a75 exhibit sudden 
changes of the structure at « 200 km from a steep de¬ 
cline to plateaus. Similar situations are also found in the 
+x and —y directions of H4-Q5a75 and —x and -\-y direc¬ 
tions of MSl-Q7a75, except for pronounced low-density 
regions between disk regions and plateaus. These gaps 
are more prominent for systems with a larger neutron- 
star radius at a fixed time (i.e., 10 ms) from the onset of 
merger and eventually disappear as tidal tails fall back. 
When material spreads in a nearly axisymmetric man¬ 
ner with (/Jej ^ Stt, plateaulike profiles are observed in all 
the directions like ALF2-Q3aO. In any case, the plateaus 
change to rapidly decaying profiles at their outer edges. 

The ejecta as an unbound portion is smoothly con¬ 
nected to a bound portion in the plateau regions. When 
the ejecta mass is large, the ejecta tends to occupy a 
large fraction of plateau material, particularly along a di¬ 
rection with the fastest expansion. The highest-density 
direction always disagrees with the fastest-expanding di¬ 
rection, in which the rest-mass density is typically lower 
by an order of magnitude at a given radius than the high¬ 
est. For example, the rest-mass density of the ejecta is 
the highest in the —y direction for APR4-Q3a75, whereas 
the fastest direction is the —x direction. This is be¬ 
cause low-density material is ejected from the outer part 
of neutron stars prior to the high-density material from 


the inner part during mass ejection driven by the tidal 
torque. The ejecta of ALF2-Q3aO is more axisymmetric 
than those of the other models, and a bump at « 650 km 
in the +y direction reflects the rear-end collision of the 
tidal tail with (^ej > 27r. Note that the spatial distri¬ 
bution of the dynamical ejecta is different from that for 
binary neutron star mergers, where a moderately steep 
power law with the index « —3.5 is observed |115] . 

The ejecta evolves in an approximately homologous 
manner. That is, the velocity of each fluid element is 
kept approximately constant, and its position and den¬ 
sity evolve according to the free-expansion law, 

r oc t , p oc (40) 


Figure [T^ shows rest-mass density and velocity profiles at 
5, 10, and 15 ms after the onset of merger in the —x and 
—y directions of APR4-Q3a75. In these plots, the radius 
and rest-mass density are scaled according to Eq. (40) 
so that those at 5 and 15 ms can be compared directly 
to those at 10 ms. Both the density and velocity profiles 
overlap approximately among different time slices after 
the scaling, and the agreement is particularly good be¬ 
tween 10 and 15 ms. These facts imply that homologous 
expansion is achieved at the late phase. We also observe 
approximate homologous expansion for other models, but 
the deviation is slightly more severe for a larger value of 
Q at a fixed time (i.e., 10 ms) due probably to stronger 
residual gravitational binding. 


2. Polar direction 

Figure [15] shows rest-mass density and velocity profiles 
along the z axis for H4-Q3a5. Because our purpose is 
to study the formation of an envelope, profiles at several 
time slices are shown together without scalings. At 5 ms 
after the onset of merger, no unbound material is found, 
and the rest-mass density is very low everywhere. This is 
because the tidal torque does not eject material toward 
the polar region. Material is pushed significantly toward 
the polar region only after the shock heating in the disk 
region sets in. This is reflected in the increase of the rest- 
mass density for t — tmerge ^ 10 ms. Unbound material 
is ejected from the disk with v « 0.3c in the beginning 
and is beyond a radius of 1000 km by « 35 ms for this 
particular model. 

A long-lived envelope is formed following the shock- 
driven disk outflow. The velocity of envelope material 
is smaller than the typical ejecta velocity, and in par¬ 
ticular, the radial velocity of bound material falls below 
0.1c at 55 ms. This suggests that the envelope is in an 
approximately stationary state at this time. Indeed, the 
rest-mass density profiles do not change very much from 
25 to 55 ms. The profile may be approximated by a power 
law, p oc with its index penv ~ 2-3. The magni¬ 

tude of the rest-mass density implies that the total mass 
of the envelope formed after the merger of H4-Q3a5 is 
much smaller than that formed after binary neutron star 
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FIG. 13. Rest-mass density in the logarithmic scale along the x and y axes at 10 ms after the onset of merger for the selected 
models. Positive and negative directions are plotted separately. Solid and dashed portions of each curve denote unbound and 
bound material, respectively. Corresponding plots in Fig. are the third row and first column for the top left (APR4-Q3a75), 
first row and second column for the top right (ALF2-Q3a(iy, fourth row and third column for the bottom left (H4-Q5a75), and 
fifth row and fourth column for the bottom right (MSl-Q7a75). 


mergers mw- This could be advantageous for a hy¬ 
pothetical jet to overcome a baryon loading problem, but 
it will not be easy to obtain a collimated jet in the ab¬ 
sence of a heavy envelope. Firm conclusions to the jet 
propagation require an extensive study of disk winds. 


It takes a long time for the remnant of a high-mass- 
ratio binary merger to develop a long-lived envelope in 
the polar region. Figure [T6| shows rest-mass density and 
velocity profiles along the z axis for MSl-Q7a75. In 
this model, the ejecta generated by the disk are beyond 
1000 km only for t — tmerge ^ 45 ms, and material behind 
it exhibits significantly more time variability than that 
of H4-Q3a5. The velocity profile with v > 0.1c also indi¬ 
cates significant time variability. It can, however, still be 
seen that the rest-mass density of the envelope is compa¬ 
rable to that of H4-Q3a5 (Fig. 15). Thus, we may safely 
conclude that the mass of the envelope formed after the 
merger is much smaller for black hole-neutron star bina¬ 
ries than for binary neutron stars unless (or possibly even 
if) binary parameters are extreme as far as the dynamical 
processes are concerned. 


3. Velocity distribution 


Figure [T^ shows the velocity distributions of dynami¬ 
cal ejecta normalized by the ejecta mass, Mej, measured 
at 10 ms after the onset of merger for selected mod¬ 
els. Namely, integrating each distribution over the ve¬ 
locity returns unity. They are derived by analyzing un¬ 
bound material on the equatorial plane as described in 
Sec. II El and we checked that estimation at different 
time slices gives very similar results. 

All the models exhibit a relatively flat distribution with 
a cutoff at low and high velocities rather than, say, a 
power-law distribution. This agrees semiquantitatively 
with previous results obtained in Newtonian simulations 
[35) . This distribution implies that the density struc¬ 
ture of ejecta can be approximated by p oc v~'^ oc 
within the range between lower and higher cutoff ve¬ 
locities, because the free-expansion law, Eq. (40), gives 
dM/dv oc pv^. This observation is largely consistent with 
the spatial profile shown in Eig. 

The velocity distribution is shifted toward larger ve¬ 
locities when the ejecta mass is larger (see the top panel 
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FIG. 14. Rest-mass density on a logarithmic scale (top) and 
radial velocity (bottom) along the —x and —y directions at 
5, 10, and 15 ms after the onset of merger for APR4-Q3a75. 
Solid and dashed portions of each curve denote unbound and 
bound material, respectively. Assnming homologous expan¬ 
sion, the radius is multiplied by 2 for the Sms profile and 
divided by 1.5 for the 15 ms profile, so that they can be com¬ 
pared directly with that at 10 ms. Similarly, the density is 
divided by 2® for the Sms profile and multiplied by (1.5)® for 
the 15 ms profile. The velocity is not scaled and is truncated 
at p = 10®gcm“® to avoid showing artificial atmospheres 
accumulated near ejecta surfaces. Truncation of profiles at 
~ 850 km for 15 ms is due to escape of material from the sec¬ 
ond largest domain, outside which hydrodynamic evolution 
equations are not solved. 


of Fig. [8 for visual comparisons). We also find that 
the distribution tends to be shifted toward larger veloc¬ 
ities when the mass ratio, Q, is larger. This is con¬ 
sistent with the observations of Mej/Mbd and Uave in 
Sec. Ill C 2[ where the dynamical ejecta from a higher- 
mass-ratio binary is seen to be energetic. Previous 
numerical-relativity simulations also found this tendency 


FIG. 15. Rest-mass density in the logarithmic scale (top) 
and radial velocity (bottom) along the z axis at several time 
slices for H4-Q3a5. Distances are also shown in the loga¬ 
rithmic scale. The time is given as that after the onset of 
merger, tmerge. Solid and dashed portions of each curve de¬ 
note unbound and bound material, respectively. The velocity 
of material with p < 10® gcm~® is not shown. 


E. Fallback 


The fallback rate as a function of time is found to obey 
a power law with the index —5/3 irrespective of the mod¬ 


els. Figure 18 shows fallback rates determined by the 
method described in Sec. |II analyzing matter profiles 
at 10 ms after the onset of merger for the selected mod¬ 
els. Aside from statistical fluctuations due to the limited 
number of grid data, overall behavior is consistent with 
the structureless power law, Mfb oc and no sig¬ 

nificant time evolution is found when we compute Mfb 
at different time slices. This power-law fallback rate is 
known to be achieved after the tidal disruption of main 
sequence stars by supermassive black holes [II611II7]. 
The same power law is found for black hole-neutron star 
binaries in Newtonian simulations [53 [50] and is also 
reported in a numerical-relativity simulation for a single 
binary model with the F = 2 polytrope |118j . Our results 
confirm their findings for a wide range of binary parame¬ 
ters in numerical relativity. Nuclear interaction neglected 
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FIG. 16. The same as Fig. but for MSl-Q7a75. 


in this study may not be important, because Newtonian 
studies show that nucleosynthesis in the nuclear statisti¬ 
cal equilibrium does not modify the power-law behavior 
[311 [no] and r-process heating can modify it only on rare 
occasions m- 

This power-law behavior implies that the mass spec¬ 
trum with respect to specific energy takes a constant pro¬ 
file, i.e., dMfh/dE = const. The usual reasoning be¬ 
hind the power-law index —5/3 is the combination of 
dMfh/dPfi:, = {dMf[y/dE){dE/dPf[y), the Keplerian rela¬ 
tion Pfb oc oc (—and the assumption that 
dMfb/dE is constant. The first and second relations are 
universal. The third assumption is verified for the tidal 
disruption of main sequence stars by various hydrody¬ 
namic simulations (e.g.. Ref. [noD and is pointed out 
to be more appropriate for a stiffer polytrope due to 
stronger shock interaction |121j . Because the neutron- 
star self-gravity cannot be neglected and shocks do not 
appear to play a significant role in energy redistribution 
for a neutron star disrupted by a stellar-mass black hole, 
the reason for constant energy distribution is nontrivial 
and may be worth future investigation. 

Although the overall magnitude of the power law is not 
computed very accurately by our approximate estimation 
method, we may safely conclude that the fallback rates 



FIG. 17. Velocity distribution normalized by the ejecta mass 
measured at 10 ms after the onset of merger for selected 
models. The velocity is dehned as -y/l — l/(—ut)^. We use 
dMej/dv rather than its integration over a finite velocity in¬ 
terval, because the former quantity is independent of binning. 



t(s) 

FIG. 18. Fallback rate measured by analyzing matter profiles 
at 10 ms after the onset of merger for selected models. A 
power law 10 “^Mq s“^( t/1 is also plotted (black dot- 

dashed line) as an eye guide. Apparent deviation from the 
power law at t > 1 s is ascribed to the limited number of grid 
data, and the power law is recovered if we compute Mfb using 
a wide time interval. 

span Mfb 10“"^-10 “^Mq s“^(t/l when substan¬ 

tial mass ejection occurs. Because the periapsis distance 
of the fallback material is found to agree approximately 
with the radius at which the neutron star is disrupted, the 
material will join the accretion disk before reaching the 
periapsis. Thus, the black-hole accretion rate and elec¬ 
tromagnetic luminosity could be smaller than the fallback 
rate (see Refs. [OTlIlOb] for relevant discussions). 

In this analysis, the center of mass is always assumed 
to be located at the coordinate origin. This is not jus¬ 
tified in a rigorous manner, because the remnant black 
hole-disk system acquires a substantial velocity of 0(100) 
kms“^ by two mechanisms. One is backreaction from the 
anisotropic mass ejection [331111, and the other is recoil 
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due to the anisotropic gravitational-wave emission m- 
We will describe the former and latter in Secs. ImFl and 
|IIIG[ respectively. 


F. Remnant disk and black hole 


Because remnant disks and black holes are thoroughly 
investigated in previous work [5S] , we describe their prop¬ 
erties only briefly. The amount of mass outside the ap¬ 


parent horizon, Mr^rAH i is shown in Table III and is 
discussed in Sec. |III C 2[ Typical accretion time scales 
due to purely hydrodynamic processes are estimated to 
be 30-300 ms when measured at « 10 ms after the onset 
of merger irrespective of the models. We do not go into 
details of accretion dynamics, expecting that realistic be¬ 
havior will be determined by unincorporated physics like 
neutrino processes and magnetohydrodynamics. 

One feature of remnant disks overlooked in our previ¬ 
ous studies is the existence of standing spiral accretion 
shocks. Figurej^shows rest-mass density profiles on the 
equatorial plane in the central region at different time 
slices for H4-Q5a75. This figure (see also the right panel 
of Fig. shows that sharp spirals extending to the appar¬ 
ent horizon stay in approximately the same location over 
10 ms without exhibiting any rotation. Similar structures 
are found for most of the models as long as the remnant 
disk is appreciable, and we find no cases in which this spi¬ 
ral structure disappears by the end of simulations, which 
is at > 50 ms after the onset of merger for the longest 
runs. The standing spiral shocks appear to be formed as 
a trace of self-collision of tidal tails rather than as a re¬ 
sult of disk instability. This spiral structure should serve 
to dissipate the angular momentum of disk material and 
enhance mass accretion by the remnant black hole. 

The mass and dimensionless spin parameter of the rem¬ 
nant black holes at 10 ms after the onset of merger are 
listed in Table [TV] They are consistent with our previous 
results for models with comparable binary parameters 
15^ 155]. After the measurement, the dimensionless 
spin parameters increase by up to « 0.01 due to long¬ 
term accretion depending on the models. Thus, the val¬ 
ues of xi shown in Table [IV| may be regarded as the lower 
limits of hypothetical final configurations which will be 
achieved by purely hydrodynamic processes. 

The remnant black hole-disk system including fallback 
material receives a recoil velocity due to the backreac- 
tion of anisotropic mass ejection gaiig, which we call 
the ejecta kick velocity. The ejecta kick velocity, 14j, is 
estimated by linear-momentum conservation as 




rei ~ 


.21 

Too 


= 555 kms 




O.OIMq 


Q + 1 


(41) 


where Mns = 1 .35M0 is assumed. For simplicity, the 
mass of the remnant black hole-disk system is approxi¬ 
mated by mo in this expression, neglecting energy loss 


to the ejecta and to gravitational waves. The former is 
< 0.02mo, and the latter is < 0.03mo for the cases con¬ 
sidered here, where the energy radiated during the very 
early inspiral phase that existed before our initial condi¬ 
tion, mo — Mo, is also taken into account. Because the 
ejecta mass is large only when tidal disruption occurs 
at a distant orbit and gravitational radiation is not very 
strong, the sum of both does not exceed 0.03mo. 

Values of the ejecta kick velocity for each model are 
presented in Table IV This table shows that t4j can 
be several hundreds of km s“^ when mass ejection is ef¬ 
ficient and easily dominates kick velocities due to the 
gravitational radiation reaction, Vgwj which we discuss 
in Sec. IIIIGI 


G. Gravitational waves 


Gravitational waves from black hole-neutron star bi¬ 
naries are thoroughly investigated in our previous work 
15^ 135] , and derived waveforms are used to construct 
phenomenological models aiming at data analysis pm 
1125] . In the following, we instead discuss integrated or 
instantaneous properties of gravitational waves. 

The energy, linear momentum, and angular momen¬ 
tum carried away by gravitational waves are presented 
in Table |IV| While the energy AEqw and the angular 
momentum A Jqw are presented as they are, the magni¬ 
tude of linear momentum APqw is shown instead as the 
velocity imparted to the remnant black hole-disk system 
including fallback material. 


Vgw 


APqw 

Too 


(42) 


where we adopt mo as in Eq. (41). We call Vow the 
gravitational-wave kick velocity. Although the accuracy 
in computing APqw is not very high due to mode cou¬ 
plings, we do not find Vgw larger than lOOkms”^ for 
the models considered in this study. Broadly speaking, 
the ejecta kick velocity, t4j, dominates the gravitational- 
wave kick velocity, Vgwj when Mej > O.OIMq. 

The ejecta kick velocity and gravitational-wave kick 
velocity partially cancel out each other, because their 
angles <i)ej and $gw point in approximately opposite di¬ 
rections. Figure pOj shows the difference between them, 
$ej — 4)GW, vs Mej. The differences cluster around tt ir¬ 
respective of the model parameters, and this means that 
ejecta and gravitational waves carry linear momenta in 
opposite directions. This tendency does not depend on 
grid resolutions. While the origin of anticorrelation is 
nontrivial, it is reasonable that 4>ej — $gw prefers a spe¬ 
cific value, because both dynamical mass ejection and 
linear-momentum emission are determined primarily by 
merger dynamics including tidal disruption. The largest 
velocity in the coalescence event is achieved by the plunge 
motion of material promptly swallowed by the black hole 
after the tidal disruption, and the plunge should emit 
the linear momentum efficiently in its direction due to 
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FIG. 19. The same as Fig. (H4-Q5a75) but at different time slices. 


TABLE IV. Characteristic physical quantities associated with the remnant black hole measnred at 10 ms after the merger 
and with gravitational waves for our fiducial, V = 60 runs. Mbh,! and xi the mass and dimensionless spin parameter, 
respectively, of the remnant black hole. Vej and Vbw are the magnitude of velocities imparted to the remnant black hole-disk 
system due to the ejecta backreaction and gravitational-wave recoil, respectively. AEaw and A Jgw are the energy and angular 
momentum, respectively, radiated by gravitational waves. 


Model 

MuH.f [Mq] 

Xi 

Vej (kms 

Vgw (kms 

AEgw[Mq] 

AJgw[-A^0] 

APR4-Q3a75 

5.07 

0.87 

100 

90 

9.3 X lO-"" 

5.6 

ALF2-Q3a75 

5.02 

0.86 

500 

60 

6.3 X 10-2 

4.6 

H4-Q3a75 

4.99 

0.88 

500 

60 

4.9 X 10-2 

4.0 

MSl-Q3a75 

4.97 

0.88 

800 

20 

4.0 X 10-2 

3.5 

APR4-Q3a5 

5.17 

0.77 

20 

70 

9.2 X 10-2 

5.3 

ALF2-Q3a5 

5.10 

0.76 

300 

70 

6.1 X 10-2 

4.3 

H4-Q3a5 

5.07 

0.76 

300 

50 

4.8 X 10-2 

3.7 

MSl-Q3a5 

5.05 

0.75 

600 

50 

3.9 X 10-2 

3.3 

APR4-Q3aO 

5.26 

0.55 

< 1 

60 

8.2 X 10-2 

4.3 

ALF2-Q3aO 

5.26 

0.56 

20 

30 

6.1 X 10-2 

3.8 

H4-Q3a0 

5.20 

0.55 

70 

40 

4.6 X 10-2 

3.2 

MSl-Q3aO 

5.16 

0.53 

200 

40 

3.6 X 10-2 

2.8 

APR4-Q5a75 

7.80 

0.85 

30 

20 

0.16 

10 

ALF2-Q5a75 

7.69 

0.83 

400 

40 

0.11 

9.1 

H4-Q5a75 

7.65 

0.83 

400 

70 

9.0 X 10-2 

8.0 

MSl-Q5a75 

7.62 

0.83 

700 

50 

7.5 X 10-2 

7.2 

APR4-Q5a5 

7.90 

0.71 

< 1 

30 

0.13 
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the large velocity (see Refs. |12611127] for relevant dis¬ 
cussions). A possible explanation of the anticorrelation 
between $ej and 4 >gw is that the linear momentum is 
emitted right after the tidal disruption primarily in the 
direction of the plunge motion, which should be oppo¬ 
site to the ejecta motion. This anticorrelation implies 
that the realistic value of the remnant velocity is given 


approximately by |I4j — VqwI- 

Finally, we comment on the possible existence of a 
tight correlation between the strength of tidal effects 
and gravitational-wave frequency at the maximum ampli¬ 
tude, which is suggested to exist for binary neutron stars 
[12811129] . Figureshows a dimensionless gravitational- 
wave frequency of the (2, 2) mode at the maximum am- 
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FIG. 20. Difference between the angle of ejecta linear mo¬ 
mentum, i&ej, and of gravitational-wave linear momentum, 
3>gw, vs the ejecta mass, Mej. We restrict the range to 
Mej > O.OOIM©. 
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FIG. 21. Dimensionless gravitational-wave frequency of (2, 2) 
mode gravitational waves, moIl 22 ,peak, at the amplitude peak 
as a function of the tidal coupling constant k defined by 
Eq. (43l. The dashed line is a fit obtained from binary neu¬ 
tron star simulations due to Ref. 


plitude, TOo^ 22 ,peak) as a function of a tidal coupling con¬ 
stant, 


2Q k _ 3Q 

(i + Q)5^ “ OTW ’ 


(43) 


adapted to black hole-neutron star binaries |130] (see also 
Ref. |131j b This figure suggests the existence of relations 
independent of the mass ratio and equation of state. If 
the correlations are tight, it implies that the finite-size 
effect in the black hole-neutron star binary merger is de¬ 
scribed fairly well by the quadrupolar tidal deformability 
up to tidal disruption. These relations depend on the 
black-hole spin, and our results suggest that mof^ 22 ,peak 
is smaller for a larger black-hole spin. This agrees with 
Ref. m- The same value of x rnay not be compared 
directly among different mass ratios except for nonspin¬ 


ning cases, and effective spin parameters weighted by the 
mass ratio such as x[l + 3/(4(5)]Q^/(l + (see, e.g.. 
Refs. |132l 1133] 1 will be more appropriate. If such cor¬ 
relations are confirmed accurately by future simulations, 
they would help to extract neutron-star equations of state 
without detailed analysis of the phase evolution, just 
as knowledge of the cutoff frequency would do [H] [5D1 - 

mm- 

Because the gravitational-wave amplitude peaks dur¬ 
ing the rapid increase of the frequency, the error of 
™o^ 22 ,peak is not very small. Typical errors are esti¬ 
mated to be « 5% due to eccentricities« 10% due to 
the finite resolution, and « 5% due to the gravitational- 
wave extraction method like extraction radii (even with 
the extrapolation). Hence, the total error may be « 20% 
in the worst case. 

A relation satisfied by nonspinning black hole-neutron 
star binaries (if really exists) is not necessarily the same 
as that by binary neutron stars, because the merger 
dynamics is very different. We include a fitting curve 
derived from binary neutron star simulations m in 
Fig. We cannot determine whether the relations are 
different or not from the current data by two reasons. 
One is the numerical error associated with each simu¬ 
lation. The other is the fact that tidal coupling con¬ 
stants, K, spanned by binary neutron star simulations are 
much larger than those by black hole-neutron star binary 
simulations, and thus the extrapolated relations cannot 
be seriously trusted. Specifically, the relation derived 
in Ref. |128j is obtained by fitting results of simulations 
with 26 < K < 440, none of which overlaps with that 
in our current simulations. It may be worth future in¬ 
vestigation to test whether relations are distinct between 
binary neutron stars and nonspinning black hole-neutron 
star binaries. 


IV. ELECTROMAGNETIC COUNTERPART 

In this section, we discuss expected characteristics of 
electromagnetic counterparts based on the properties of 
dynamical ejecta derived by our simulations. We focus 
primarily on the effect of anisotropy, which is character¬ 
ized by the opening angle in the equatorial plane, ipej, and 
in the meridional plane, 9ej, on the macronova/kilonova 
and synchrotron radio emission [231 US] ■ A. con¬ 
cise summary of the main results derived in this section 
is found in Ref. (35) j in which other aspects of the ejecta 
like gravitational-wave memory emission and cosmic-ray 
acceleration are also discussed. 

For simplicity, we adopt slightly different notations for 
ejecta quantities in this section from those in other sec¬ 
tions. Specifically, we denote the ejecta mass by M in- 


This is estimated as twice the eccentricity in the inspiral phase, 
because it is difficult to isolate the eccentricity contribution dur¬ 
ing the merger phase. 
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stead of Mej. The opening angles are denoted by 9 and (p 
instead of 0ej and ipej, respectively. Recall that 6{= dej) 
is defined as the half-opening angle taking the equatorial 
symmetry into account, and a full sphere corresponds to 
0 = 7r/2 and p = 2 tt. We also adopt short-hand nota¬ 
tions M _2 = M/(O.O1M0), 6i = 6/{l/5), and pi = pj'K. 
We recover the speed of light, c, everywhere. 


A. Macronova/kilonova 

The macronova/kilonova is quasithermal radiation 
from the ejecta heated by the decay of unstable r-process 
elements. The dynamical ejecta from mergers of black 
hole-neutron star binaries will be composed primarily of 
neutrons as discussed in Sec. |IIIB 1[ and then r-process 
elements should be synthesized [5S1 [27j . After the neu¬ 
trons are exhausted within a few seconds, /3-decay and fis¬ 
sion of unstable r-process elements heat the ejecta[^ The 
heated ejecta emit radiation primarily in red-optical and 
infrared bands on a day-to-month time scale |135) , where 
a bunch of Doppler-broadened lines associated with the 
complicated energy-level structure of r-process elements 
blanket the emission in blue-optical and ultraviolet bands 


1. Analytic model 


Qualitative features of the macronova/kilonova from 
the anisotropic ejecta can be understood by modifying 
the prototypical model for spherical ejecta proposed in 
Ref. [20]. In this section, we introduce short-hand no¬ 
tations V-i = R/(0.1c) for the surface velocity, ki = 
K/(10g“^ cm^) for the opacity, and /_6 = //10“® for the 
heating efficiency. The precise meaning of these quanti¬ 
ties is explained in the following. 

We approximate the hydrodynamic evolution of the 
ejecta by the free expansion of a uniform-density trun¬ 
cated sphere characterized by the opening angles 6 and 
p. The radius of the (truncated) sphere is given by 
R{t) = Vt using the surface velocity V, and thus the 
rest-mass density of the ejecta is 


SttM 

^ A9pVH^ ■ 


(44) 


In this uniform-density free-expansion model, the surface 
velocity is related to the aver age v elocity of the ejecta de¬ 
fined by Eq. (20) via V = i/S/Suave ~ l-3uave- Because 
the ejecta material is expected to be radiation-dominated 
in the relevant epoch due to r-process heating, the inter¬ 
nal energy density u is related to the pressure P and 


Some of the energy liberated in the /3-decay does not contribute 
to the heating because of the energy deposited to neutrinos and 
7 -ray photons m- The latter does not escape freely in the early 
stage of the ejecta evolution and contribute to the ejecta heatup. 


temperature T by w = 3P = where a is the radia¬ 
tion constant. The time evolution of the internal energy 
density is derived by the first law of thermodynamics as 


r,du 2 37rM . 37r 

^ ~di^ ^ ^ AOpV^^ ~ ABpV^^' 


(45) 


where e is the specific heating rate and L is the luminos¬ 
ity. Time-dependent quantities are tt, e, and L. 

We assume that the specific heating rate is given by a 
power law e(t) = fc^/t parametrized by the heating effi¬ 
ciency / in the same manner as the spherical model [20j . 
An appropriate value of the heating efficiency, /, will de¬ 
pend significantly on the electron fraction (the number 
of electrons per baryon ) unnidns]. The uncertainty is 
particularly high when fission is an important heating 
source rather than /3-decay of elements near the stability 
line |100) . In this study, we take the fiducial value of / 
to be 10“® following Ref. |l()()j . 

We give the luminosity by a diffusion approximation 
in a similar manner to the spherical model m but as¬ 
suming geometry adapted to anisotropic mass ejection. 
The assumption is that the radiation is emitted not from 
the truncated spherical surface but from the cross section 
of truncation. In the language of our simulations, pho¬ 
tons from the anisotropic ejecta are assumed to escape 
mainly into the ±z directions, and the emitting surfaces 
are taken to be those observed from the ±z direction 
like ones depicted in Fig. [^ The temperature gradient 
dT/dr relevant to the diffusion flux is approximated by 
« T/{9R) rather than « T/R of the spherical ejecta un¬ 
der this assumption. Thus, the flux may be given by 


F{t) 


<JsbT* 
Kp6R ’ 


(46) 


where k is the opacity and ctsb is the Stefan-Boltzmann 
constant. In this estimation, a factor of order unity is 
neglected in exactly the same manner as in Ref. [20] ■ The 
emitting area is then given approximately by 2xpR'^(2 = 
pR^, where the first “2” stands for two emitting surfaces 
at +z and —z, and therefore the bolometric luminosity 
may be given by 


m 


3^' " 


(47) 


This expression does not reduce to that for the spherical 
ejecta even if we adopt 0 = 7r/2 and p = 2 tt because of 
different assumptions. The neglected truncated spherical 
surface has the area {A/tt)9pR^^ and thus the luminosity 
may be underestimated by a fraction of (4/7r)0 « 30%. 
Although this term can be included with no difficulty, we 
omit this contribution so that the parameter dependence 
becomes clear. 

The value of opacity, k, is highly uncertain due to our 
incomplete knowledge of r-process elements and their 
line features [iMllm] . Although the realistic opacity 
of r-process elements is safely assumed to be dominated 
by various bound-bound transition lines in optical and 
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ultraviolet wavelengths, no complete line list exists so 
far. In this study, we take the fiducial value of k to 
be 10g“^cm^, because this approximately reproduces 
results obtained by radiation transfer simulations per¬ 
formed adopting currently available line lists |137j . The 
gray approximation adopted in this model is not realistic 
and limits the predictability of the spectra. 

The thermodynamic evolution equation can be solved 
analytically. For this purpose, it is convenient to cast the 
equation into a dimensionless form. First, we normalize 
the surface velocity by the speed of light as /3s = Vjc 
in the usual manner. Next, we define a characteristic 
time scale by the condition that the optical depth of a 
characteristic path becomes unity, npOR = 1, and this 
gives a critical time of the onset of transparency. 


tc 


SttkM 


4ipV^ ■ 


(48) 


Finally, a characteristic internal energy density can be 
defined by 


Introducing dimensionless variables t = tjtc and u = 
m/uc, we obtain the dimensionless evolution equation. 


du /4 Siri \ _ 1 

di \i 166*/3s ) i'^ ’ 

which has an analytic solution 


(50) 


u{t) = ^ exp 


SttP \ l32ePs ly( I Stt 
329Sir W 320/3sJ ’ 


(51) 

where C is the integration constant and Y is Dawson’s 
integral defined by 


Y{x) 



(52) 


Because the initial internal ener gy o f the ejecta is rapidly 
lost due to the adiabatic coolingwe may safely set the 
integration constant, C, to be zero. The key issue which 
allows us to derive this analytic solution for nonspheri- 
cal ejecta is that the temporal dependence of each term 
(adiabatic cooling, radiative cooling, and heating) is not 
affected by the geometry in our model. 

The peak time, peak bolometric luminosity, and ef¬ 
fective temperature at the peak time can be estimated 
using this solution. Note that Dawson’s integral takes 


the maximum value yp « 0.54 at Xp « 0.92. The peak 
time is 


^peak — 


I8k0M 


(fVc 
= 11 day 




(53) 


The peak bolometric luminosity is 


Tpeak — Vp 
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= 1.8 X 10^°ergcm-3 . 

(54) 

The effective temperature is defined from the diffusion 
flux, Eq. (461, by Teff = (F/ctsb)^^"^, and its value at the 
peak time is 


Tpeak 


vl'U 

\8a?9^K^Mv) 

= 1900 K (55) 


These expressions share the same parameter dependence 
as those derived in Ref. [49j by random-walk arguments. 
This indicates that the parameter dependence is robust 
as far as similar assumptions are adopted. 

In the range of opening angles observed in our nu¬ 
merical simulations, the macronova/kilonova from black 
hole-neutron star binaries tends to peak slightly earlier 
with slightly higher bolometric luminosity than that from 
the spherical ejecta for given values of other parameters. 
This tendency is also observed in radiation transfer sim¬ 
ulations |135j . When the opening angle in the meridional 
plane, 9, is small, the peak time becomes early and the 
peak bolometric luminosity increases. The reason for this 
is that photons can escape easily from the ejecta when 9 is 
small. Specifically, the optical depth Kp9R at tpeak is pro¬ 
portional to 9~^ and independent of ip. When the open¬ 
ing angle in the equatorial plane, ip, is small, the peak 
time becomes late and the peak bolometric luminosity 
decreases. The reason is that a small value of p increases 
the rest-mass density, optical depth, and characteristic 
time scales. The dependence of luminosity may be un¬ 
derstood by the constancy of Tpeak^peak for both cases. 
The combined effects of these two angles tend to prefer 
the slightly earlier peak with slightly brighter emission. 
When the ejecta becomes transparent, the bolometric lu¬ 
minosity does not depend on the geometry, because we 
simply have L = eM even within this model derived with 
the diffusion approximation!^ 

At a given time, the material temperature T and effec¬ 
tive temperature Teff are higher for the anisotropic ejecta 


Further energy injection could modify the thermodynamic evo¬ 
lution via different e{t) m- 


It is probable that the photons are depleted when the free-free 
emission becomes inefficient and this effect is not taken into 
account in the current model. 
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than for the spherical one due to different geometry. In 
typical situations, the material temperature, T, is higher 
by about a factor of 2, and this agrees approximately 
with the result of Ref. |135j . The reason of the high tem¬ 
perature is that the decay heat of unstable r-process ele¬ 
ments is deposited to a small volume for a given mass and 
velocity of the ejecta. Accordingly, the effective temper¬ 
ature, Teff, is higher by « 30%-50%. Even this amount 
of difference could have a significant effect on the ob¬ 
served flux [not to be confused with the diffusion flux, 
Eq. (46), which is trivially related to Teg] in optical and 
near-infrared bands, because the typical value of Tpeak is 
in the infrared band. Thus, a small increase of Teg could 
enhance the flux at optical and near-infrared bands. In 
fact, we see that absolute/apparent magnitudes increase 
by 1-2 in these bands if we assume a perfect blackbody 
spectrum. We do not, however, regard this amount of in¬ 
crease as very quantitative, because realistic spectra will 
be very different from the blackbody. The high temper¬ 
ature could also affect possible dust formation |138j . 

The dependence of the peak quantities on /, k, M, and 
V is the same as that of the spherical model [2D] . A large 
value of / increases the luminosity, a large value of k de¬ 
lays the peak and decreases the luminosity, a large value 
of M delays the peak and increases the luminosity, and 
a large value of V hastens the peak and increases the lu¬ 
minosity. Among these parameters, the ejecta mass can 
become much larger for black hole-neutron star binaries 
than for binary neutron stars |44] . and higher luminos¬ 
ity could be achieved |112j . In fact, this difference may 
dominate corrections due to the opening angles. Typi¬ 
cal velocities of the ejecta cannot be very different. The 
heating efficiency and opacity could change reflecting dif¬ 
ferent compositions of the ejecta, but we do not discuss 
them in this study. 


2. Directional dependence, line, and polarization 

We briefly discuss possible aspects of the 
macronova/kilonova from anisotropic ejecta that cannot 
be captured in the analytic model developed above. An 
obvious outcome of the anisotropy is the directional 
dependence (see Refs. [1391 I140j for binary neutron 
stars). Emission should be brighter when observed from 
the direction perpendicular to the equatorial plane than 
in the equatorial plane. Specifically, the flux will be 
larger by 1/d « 3-5 near tpeak for the former situation. 
Accordingly, the light curves will exhibit different 
evolution near tpeak and become indistinguishable after 
the entire ejecta becomes transparent. These behaviors 
are observed in radiation transfer simulations m- 
Followup observations of electromagnetic counterparts 
will benefit from this directional dependence, because 
gravitational waves are emitted most strongly in the 
direction perpendicular to the equatorial plane, in which 
the macronova/kilonova will be the brightest. That is, 
observed binaries will be biased toward the brightest 


direction of the macronova/kilonova. 

A chance to observe spectral lines associated with r- 
process elements will be better for black hole-neutron 
star binaries than for binary neutron stars. In any case, 
it will be very challenging to observe such lines from 
the macronova/kilonova, because a bunch of lines are ex¬ 
pected to be significantly blended due to Doppler broad¬ 
ening in the ejecta with a large surface velocity. This 
broadening may be mitigated in the direction perpen¬ 
dicular to the equatorial plane, because the expansion 
velocity is smaller by a factor of 6 than in other direc¬ 
tions and spherical cases. Furthermore, the emission is 
expected to be the brightest in this direction. Thus, the 
macronova/kilonova associated with black hole-neutron 
star binaries would deserve detailed spectroscopic obser¬ 
vations to seek a (serendipitous) strong and isolated line 
(see also Ref. |136| for relevant discussions). 

Potential diagnostics of the anisotropic geometry is 
polarization induced by electron scattering, but the 
polarization degree is not likely to be high for the 
macronova/kilonova. If the optical depth to electron 
scattering is sufficiently high and lines do not contribute 
to depolarization significantly, the linear polarization ob¬ 
served from the equatorial plane would be 4%-5% be¬ 
cause of the highly deformed photosphere |141j . However, 
the number density of free electrons will be much smaller 
in the ejecta composed of r-process elements than in, e.g., 
the supernova ejecta near the peak luminosity. While 
the r-process elements have the mass number > 100, the 
ionization degree will not be particularly high around 
the peak of the macronova/kilonova [13611137] . Hence, 
the opacity for the electron scattering will be lower by 
about 3 orders of magnitude than that for bound-bound 
transitions if k = 10 g“^ cm^ is an appropriate represen¬ 
tative. The optical depth to electron scattering will be 
only 0(10“^) near tpeak when the total optical depth is 
at most 0(10), and therefore the polarization degree may 
be reduced by a similar factor. In addition, interaction 
with lines will further depolarize the radiation [142j . 


B. Synchrotron radio emission 


Nonthermal radiation such as synchrotron emission is 
expected to arise from blast waves formed between the 
ejecta and ambient interstellar medium in a similar man¬ 
ner to the supernova remnant and gamma-ray burst af¬ 
terglow PH [DS] . Subrelativistic blast waves will develop 
as the ejecta sweeps the interstellar medium, and the 
kinetic energy of the ejecta is converted to postshock in¬ 
ternal energy. A fraction of the internal energy at the 
forward shock will be converted to energy of nonthermal 
electrons assembled from the interstellar medium and of 
amplified magnetic fields. The accelerated electrons will 
radiate synchrotron emission in a magnetized environ¬ 
ment, and the emission would be observed in radio bands 
[ 3121 ] and possibly in optical, x-ray, and y-ray bands [25] . 
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1. Radiation at the deceleration time 

The most luminous emission is expected when the 
ejecta begins to be decelerated signihcantly, and the de¬ 
celeration time tdec depends on the ejecta geometry for 
a given mass and velocity of the ejecta. We describe the 
synchrotron radio emission expected at tdec adopting a 
simplified version of the nonrelativistic model developed 
in Ref. [55] (see also Ref. |143j b We do not, however, 
attempt to model the time evolution in this study, be¬ 
cause the lateral expansion should become important af¬ 
ter the deceleration time for the anisotropic ejecta. While 
the late-time evolution of the spherical ejecta will be de¬ 
scribed reasonably by Sedov-Taylor’s self-similar solu¬ 
tion as in Ref. [7], it is difficult to formulate the lat¬ 
eral expansion of the anisotropic ejecta in a simple man¬ 
ner. We introduce short-hand notations w_i = u/(0.1c) 
for the ejecta velocity, ng = n/(lcm“^) for the ambi¬ 
ent number density, ee,-i = ee/ 0.1 for the fraction of 
postshock internal energy given to nonthermal electrons, 
cb,-! = es/O.l for the fraction of postshock internal en¬ 
ergy given to magnetic fields, and D 2 = 13/(100Mpc) 
for the distance from the observer to the site of binary 
coalescence. We also introduce the fraction rj of accel¬ 
erated electrons and power-law index p of the Lorentz 
factor distribution. 

The ejecta is decelerated significantly when the mass 
comparable to its own is assembled from the interstel¬ 
lar medium. The deceleration radius of the anisotropic 
ejecta is given by 




= 1 - 
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(56) 


where rup is the proton mass and n is the nnmber density 
of the ambient medium. Modeling the ejecta by a single¬ 
velocity shell with v, the corresponding deceleration time 
is 


^dec — 


-^de 


= 38year Uq ■ (57) 


Here, the ejecta velocity, u, may be identihed with the 
average velocity of the ejecta, Uave, defined by Eq. (20). 


The value of the ambient density, n, should vary by orders 
of magnitude depending on the location of the binary 
coalescence, and we normalize it by the typical value of 
the Galactic disk, 1 cm“^, following previous work [ZllMl 
[5K] . These expressions reduce to values for the spherical 
ejecta when we set Q = 7r/2 and p = 27r. 

If the typical opening angles observed in numerical 
simulations are kept until the deceleration time, Rdec 
and tdec are larger by a factor of 2-3 than those for 
the spherical ejecta for given values of the other param¬ 
eters. The reason for this is that only a limited fraction 
of the volume inside i?dec is swept by the anisotropic 


ejecta. Whereas the ejecta will approach a spherical state 
to some extent before tdec [103j . the synchrotron radio 
emission from black hole-neutron star binaries will be a 
longer-lasting event than that for binary neutron stars. 

The geometry does not modify the Lorentz factor dis¬ 
tribution of nonthermal electrons and magnetic fields. 
The number of assembled electrons at t^ec is given by 

.r M 

Ne,tot — - 

TOp 

= 1.2 X 10®® M_ 2 . (58) 

Assuming that a fraction 77 < 1 of these electrons is ac¬ 
celerated to power-law distribution of the Lorentz factor 
7 e with the index p > 2 as 

(X j-P (7„ < 7e), (59) 

d'~fe 

the minimum Lorentz factor may be derived from the 
number and energy of accelerated electrons as 

Ce p - 2 rup (u/c)^ 

7m — z 

T] p — 1 TUq 2 

= 0.92 g(p)ri~^ee-iv‘ii, (60) 

where rrie is the electron mass, Cg is the fraction of the 
postshock internal energy given to the accelerated non¬ 
thermal electrons, and p(p) = (p — 2)/(p— 1). Care must 
be taken in applying this equation to subrelativistic blast 
waves, because 7 ^ can fall below unity and become un¬ 
physical, particularly when the ejecta is significantly de¬ 
celerated (not considered here). The strength of mag¬ 
netic fields is given by 

B = ^^'KeBnm-pV 

= 6.5mGauss (61) 

where es is the fraction of postshock internal energy con¬ 
verted to magnetic fields. 

Parameters characterizing microphysics, p, 77 , Cg, and 
cb, are all uncertain. Following Ref. [7], we normalize 
Eg and cb by 0.1 and take the fiducial value of p to be 
2.5, which gives g{p) = 1/3. Typical values ofp observed 
in nonrelativistic blast waves may be 2.5-3 |144j . The 
fiducial value of 77 is set to be unity. Detailed spectro¬ 
scopic observations of nonthermal radiation could deter¬ 
mine these parameters in principle [25] . 

Quantities characterizing the instantaneous spectrum 
are estimated as follows. The synchrotron frequency of 
an electron with 7 e is defined by Veile) = /(2Trmec), 
where q is the elementary charge, and the power of the 
electron is by Peije) = <^TcB‘^^^/{6Tr), where ctt is the 
Thomson cross section. The specific flux from a single 
electron at its peak frequency, Bg, is estimated to be 
Pi, « Pg/iZg = aTmeC^B/{3q) independently from the 
electron Lorentz factor. The characteristic frequency of 
the electron distribution corresponding to 7 ^ is given by 

Vm = 1.5 X 10^ Hz gip)^V~‘^el-ieB^_inl^‘^v^^. 


(62) 
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An unabsorbed specific flux, i.e., a hypothetical specific 
flux in the absence of self-absorption, at Vm is estimated 
to be 

“ 47r£|2 

= 2.4 Jy (63) 

where D is the distance from the observer to the site 
of binary coalescence. If we neglect synchrotron self¬ 
absorption and cooling, the specific flux is given by 
below Vra and above 

Vjn- We show below that the self-absorption could sup¬ 
press the radio spectrum, while the cooling is not impor¬ 
tant in the radio band. 

The self-absorption frequency Va may be obtained ap¬ 
proximately by comparing the hypothetical unabsorbed 
flux with the blackbody flux in the Rayleigh regime |143j . 
The blackbody flux at v^n is given by 

= 27rr'^7mTOe^^, (64) 

where Aem is the blackbody emitting area. While this 
should be 47ri?jgj. for the spherical ejecta, we take Aem = 
V^-^dec for the anisotropic ejecta in a similar 

manner to Sec. |IV A 1| It is readily found that Vm < 
Va when < Tb.m, and vice versa. The case that 

Vm < Va is typical for subrelativistic blast waves, and 
the self-absorption frequency is defined as the frequency 
at which the synchrotron and blackbody fluxes are equal. 
Specifically, we obtain 


Va 


f F^^rn 
V -Fi/,BB 


2/(4+p) 


= 3.9 X 10^ Hz ry-2(P-2)/(4-Hp) 

y 2 (p-1)/(4+p) (2+p)/[2(4+p)] (14-H3p)/[6(4-Hp)] 

^ ^e,-l ^B.-l ''-0 

V A.42/[3(4-hp)] (5p-2)/(4-Hp).4/[3(4+p)] -2/[3(4+p)] 

X IV1_2 V_i a- if- , 

(65) 


role, is estimated by the condition = P(7c)tdec, 

and we obtain 

bTTTOeC 

= 1.5 X 10" (66) 

The corresponding cooling frequency is given by 


Va = 4.3 X 10^2 Hz ■ (67) 


Although the cooling frequency decreases by a factor of 
several for the anisotropic ejecta due to the long deceler¬ 
ation time, this could affect the radio spectrum at high 
frequency only in a limited parameter range. 

Finally, the instantaneous spectrum for Vm < Va < Vc 
is given by 


{{ValvY)~^'^^^'^l’^{vlvYi^ {v<Vm) 

Fu _ I {ValVm)~^^~^"'^‘^{vlVaf/'^ [v^ < V < Va) 

Fu,m I (v {Va< V < Vc) 

[{Va/Vm)-^P-^^^Hv/Va)-P/^ {Va < v) 


( 68 ) 

The third segment is the most relevant to radio observa¬ 
tions, and it would be useful to reexpress the spectrum 
in this range as 


-(p-l)/2 


a = 0.12mJy(j^) 

^1-p^p-l ^(p+l)/4^(p+l)/4 


X r] 




M_2vY 


(5p-3)/2^_2 


(69) 


where the prefactor is for p = 2.5 and decreases by a 
factor of 40 as p increases from 2.1 to 30 This ex¬ 
pression indicates that the emission associated with the 
massive ejecta from black hole-neutron star binaries will 
be bright. 


2. Proper motion 


where the prefactor is given for p = 2.5 and varies 
by a factor of 2 within 2.1 < p < 3j^ The depen¬ 
dence on the opening angles is inherited from the emit¬ 
ting area, and Va is smaller by a few tens of 

percent for the anisotropic ejecta than for the spheri¬ 
cal one. The self-absorption frequency can increase to 
~ 1 GHz in a plausible parameter range, and thus the 
self-absorption could be important at low-frequency ra¬ 
dio bands for such cases. When Va < Vm, we instead 
obtain Va = {Fa^m/F^^BB)^^^Vm- 

The cooling Lorentz factor 7 c at the deceleration time, 
above which the radiative energy loss plays a significant 


Aside from the expansion, the anisotropic ejecta from 
black hole-neutron star binaries exhibits center-of-mass 
motion, and thus the proper motion of radio images 
could be observed [35]. The characteristic distance of 
the center-of-mass motion may be given approximately 
by Rcom — 4^ej^dec — Rdec(^^ej/^^ave)■ The projected dis¬ 
tance on the celestial sphere should be smaller by a factor 
of « 2 due to the angular average, whereas the observa¬ 
tional bias due to the directional dependence of grav¬ 
itational radiation should mitigate this decrease. The 
expected amount of projected travel distances is 0(1) 
pc [see Eq. (56)], and we expect the radio image of the 
anisotropic ejecta to move 0(1) milliarcsecond during its 


4® More precisely, the prefactor (6.0 X x 

l.SxlO^Hz g(p)2{p-l)/(4+p) ig applicable to all the val- The prefactor (6.5 X ^ 24 Jy g(p)P“4 jg applicable to 

ues of p > 2. all the values of p. 
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bright emission for a event at 0(100) Mpc. This amount 
of proper motion could be resolved by current radio in¬ 
struments depending on the parameters and observed fre¬ 
quency m and could help to distinguish black hole- 
neutron star binaries from binary neutron stars only by 
electromagnetic counterparts. 


V. SUMMARY 

We performed numerical-relativity simulations of black 
hole-neutron star binary mergers to study dynamical 
mass ejection. The mass ratio, black-hole spin, and 
neutron-star equation of state were systematically var¬ 
ied to reveal the dependence of the ejecta properties on 
these parameters. We found that dynamical mass ejec¬ 
tion is driven primarily by the tidal torque exerted from 
black holes to elongated neutron stars, and this process 
progresses over « 2 ms after the onset of merger. The 
dynamical ejecta is concentrated around the equatorial 
plane with a half opening angle of 10°-20° and sweeps out 
about a half of the plane, except for cases that mass ejec¬ 
tion is inefficient. Because of this anisotropy, the ejecta 
carries a bulk linear momentum, and thus the remnant 
black hole-disk system receives an ejecta kick velocity due 
to the backreaction. 

The ejecta mass can be as large as ~ O.IMq, and the 
average velocity of the ejecta defined from the kinetic en¬ 
ergy is typically 0.2-0.3c. Dynamical mass ejection tends 
to become efficient when the neutron-star compactness is 
small, the mass ratio is small, and/or the black-hole spin 
is large. The dependence of ejecta properties on the com¬ 
pactness, however, is not as simple as that of the total 
mass remaining outside the apparent horizon. This sug¬ 
gests that not only the compactness but also detailed 
properties of the equation of state influence the ejecta 
properties significantly. Furthermore, the dependence on 
the mass ratio is not always monotonic. The ratio of the 
ejecta mass to the bound mass is large when the mass ra¬ 
tio is large, and the average velocity of the ejecta is also 
large for such cases. These suggest that the dynamical 
ejecta from higher-mass-ratio binaries is more energetic 
for a given ejecta mass. 

We also found that the bound envelope along the po¬ 
lar axis of the central remnant is not as heavy as that 
for binary neutron star mergers as far as the dynamical 
processes are concerned. This would be advantageous 
for a hypothetical gamma-ray burst jet to overcome the 
baryon loading problem, while how to collimate it is un¬ 
certain in the absence of a heavy envelope. Fallback rates 
of bound material obey the canonical —5/3 power law. 
The remnant disk exhibits a standing spiral shock struc¬ 
ture, which enhances the mass accretion. 

Because the gravitational-wave kick velocity imparted 
to the remnant does not exceed 100kms“^ for our mod¬ 
els, the ejecta kick velocity dominates motion of the rem¬ 
nant. We found that ejecta and gravitational waves usu¬ 
ally carry the linear momentum in the opposite direction. 


and thus these two kick velocities would partially cancel 
out. Tight correlations between the gravitational-wave 
frequency at the maximum amplitude and tidal coupling 
constant were suggested to exist in a similar manner to 
that found for binary neutron stars. The relations for 
black hole-neutron star binaries depend on the black-hole 
spin. 

Properties of electromagnetic counterparts were dis¬ 
cussed based on the results of numerical simulations fo¬ 
cusing on the effect of ejecta anisotropy. An analytic 
model of the macronova/kilonova shows that both the 
material and effective temperatures become high for the 
anisotropic ejecta from black hole-neutron star binaries. 
We also found that the peak time is slightly early and the 
peak bolometric luminosity is slightly high for the typical 
ejecta opening angles. The synchrotron radio emission is 
long lasting for the anisotropic ejecta, and the proper 
motion of the radio images could also be observed. The 
most significant difference from electromagnetic counter¬ 
parts associated with binary neutron stars would come 
from different ejecta masses for both emission models. 
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Appendix A: Convergence and uncertainty 

Ejecta are only a fraction of the material remaining 
outside the horizon, which itself is only a fraction of a 
neutron star. Therefore, quantities associated with the 
ejecta could entail large fractional errors. Furthermore, 
our numerical simulations have various parameters both 
physical and unphysical. In this appendix, we estimate 
errors and uncertainties in our computations. We also 
discuss seemingly spurious high-velocity ejecta found in 
Sec. lHIBlI 
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FIG. 22. Merger time normalized by the total mass, 
tmerge/mo, VS grid resolutions for selected models. We as¬ 
sume a hypothetical convergence order 2.5 for all the models. 
Actual numerical data for APR4-Q3a75, MSl-Q5a75, and H4- 
Q7a75 show convergence orders 2.2, 2.9, and 3.1, respectively. 


1. Convergence with respect to the grid resolution 


Finite grid resolutions are obvious sources of errors. 
First of all, we demonstrate that reasonable convergence 
behavior is observed in our numerical simulations. Fig- 
I 22 shows the merger time, tmerge, as a function of grid 


ure 


resolutions represented by N (see Sec. IIC for the defi¬ 
nition) for selected models. The exact convergence order 
estimated from numerical data varies among models and 
typically lies between 2 and 3. Taking the different ac¬ 
curacies for different parts of our code SACRA [55] into 
account, the observed behavior is reasonable. 

Table [V] compares characteristic quantities among dif¬ 
ferent grid resolutions for selected models. It is evident 
that these quantities are not always monotonic with re¬ 
spect to the grid resolution. Such behavior is frequently 
seen in hydrodynamic quantities, which severely suffer 
low convergence order when shock waves exist. Rela¬ 
tive errors are smaller for ah than for Mgj, and this 
suggests that the accurate determination of boundaries 
separating bound and unbound material is an important 
but difficult task. If we assume the first-order conver¬ 
gence between TV = 60 and 48 results, the worst-case 
error with TV = 60 results are « 30%, 40%, and 100% 
for Mr>r-AH; -^bd) and Mej, respectively. The accuracy 
of Mej is especially low when the ejecta mass is as small 
as Mej < O.OIM 0 , and the error decreases to < 20% for 
more massive ejecta. It is reasonable that the relative er¬ 
ror is large for the small mass ejecta, where the absolute 
error is always estimated to be » 0.005-0 .OIMq for Mej. 


2. Effect of an artificial atmosphere 

An artificial atmosphere affects the ejecta properties. 
Some portion of the atmosphere happens to satisfy the 


unbound criterion, ut < —1, as a result of hydrody¬ 
namic interaction, and this error spuriously increases the 
amount of the ejecta. At the same time, the atmosphere 
decelerates physical material ejected from neutron stars, 
and this error spuriously decreases the amount of ejecta. 
Low atmospheric density will mitigate both these errors. 
Another source of error is a steep density gradient at the 
neutron-star surface, which induces spurious shock heat¬ 
ing in numerical simulations and helps the material to 
become unbound. Although this error will be suppressed 
as grid resolutions are improved to resolve the stellar sur¬ 
face accurately, lowering the atmospheric density at a 
fixed resolution does not always suppress it, because the 
shock could become strong. 


Table VI compares characteristic quantities obtained 
with different values of /at. We find that ejecta quanti¬ 
ties like Mej, Tej, and Pgj increase by « 0.5% for = 60 
when the atmospheric density is decreased by an order 
of magnitude, while the corresponding change (either in¬ 
crease or decrease) is « l%-2% for N = AO. Although 
the dominant mechanism responsible for the error is not 
certain, this suggests that the error associated with the 
artificial atmosphere will decrease significantly as grid 
resolutions are improved probably due to the suppres¬ 
sion of spurious shocks at the stellar surface. We also 
find a similar amount of variations when we change the 
value of the atmospheric power-law index riat from 3 to 
2 . 


3. Effect of thermal correction Fth 

Dynamical mass ejection is expected to be governed 
basically by zero-temperature equations of state, because 
shock heating does not play a significant role. How¬ 
ever, the ejecta properties depend weakly on the finite- 
temperature part of equations of state due to the spuri¬ 
ous shock heating. Thus, the dependence of results on 
Fth also requires investigation. 

Table |Vn] compares characteristic quantities obtained 
with different values of Fjh for selected models. Results 
do not depend monotonically on the value of Fth. This 
fact suggests that the effect of Fth is not very physical, 
and the difference is ascribed to numerical errors such 
as spurious shock heating. We checked that the differ¬ 
ences develop during dynamical mass ejection which pro¬ 
gresses over « 2 ms after the onset of merger, and late¬ 
time physical shock heating in the disk region does not 
introduce significant differences. The difference of results 
among different values of Fth is as large as « 20% when 
Mej < O.OIM 0 and tends to become small when the mass 
ejection is efficient. We regard the difference observed in 
Table [VTTI as an estimate of systematic uncertainty, which 
could converge as grid resolutions are improved. 
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TABLE V. The same as Table [III| but for runs with different grid resolutions for selected models. The fiducial grid resolution 
adopted in the body text is = 60. 


A 

^r>rAH [-^©] 

Afbd[Afo] 

MejlMe] 

Tej (erg) 

Tej [Mg] 

“^ave 

Uej 

APR4-Q3a75 

60 

0.194 

0.182 

0.0111 

5.48 X lO'’" 

2.08 X 10“® 

0.235 

0.187 

48 

0.210 

0.201 

0.0087 

3.98 X 10®° 

1.61 X 10~® 

0.227 

0.186 

40 

0.206 

0.198 

0.0079 

3.40 X 10®° 

1.41 X 10“® 

0.219 

0.178 

APR4-Q5a75 

60 

0.068 

0.059 

0.0084 

4.77 X 10®° 

8.30 X 10“^ 

0.252 

0.099 

48 

0.070 

0.063 

0.0067 

3.60 X 10®° 

7.87 X 10“^ 

0.246 

0.118 

40 

0.074 

0.069 

0.0047 

2.43 X 10®° 

6.51 X 10“"® 

0.239 

0.137 

MSl-Q5a75 

60 

0.356 

0.277 

0.0785 

5.55 X 10®® 

1.80 X 10"^ 

0.281 

0.223 

48 

0.361 

0.282 

0.0795 

5.26 X 10®® 

1.80 X 10~^ 

0.272 

0.226 

40 

0.370 

0.290 

0.0797 

5.28 X 10®® 

1.77 X 10“^ 

0.272 

0.222 

H4-Q7a75 

60 

0.194 

0.157 

0.0375 

2.78 X 10®® 

7.19 X 10“® 

0.288 

0.192 

48 

0.207 

0.171 

0.0360 

2.64 X 10®® 

6.89 X 10“® 

0.287 

0.191 

40 

0.214 

0.179 

0.0341 

2.41 X 10®® 

6.64 X 10“® 

0.281 

0.195 


TABLE VI. The same as Table [nl| but with different values of /at for APR4-Q3a75 with A = 60 and 40. The fiducial value of 
/at adopted in the body text is 10“^^. 


/at 

-^r>rAH [-^ol 

Mbd[M0] 

M.AMq] 

Tej (erg) 

Tej[M0] 

'^ave 

Vej 

O 

II 

10~®® 

0.196 

0.185 

0.0111 

5.45 X 10®° 

2.07 X 10"® 

0.235 

0.186 

10"®^ 

0.194 

0.182 

0.0111 

5.48 X 10®° 

2.08 X 10“® 

0.235 

0.187 

10"®® 

0.194 

0.183 

0.0112 

5.51 X 10®° 

2.08 X 10“® 

0.235 

0.186 

A = 40 

10-®u 

0.207 

0.199 

0.0080 

3.44 X 10®° 

1.43 X 10"® 

0.219 

0.178 

10"®® 

0.209 

0.201 

0.0078 

3.33 X 10®° 

1.38 X 10“® 

0.219 

0.177 

10"®® 

0.206 

0.198 

0.0079 

3.40 X 10®° 

1.41 X 10“® 

0.219 

0.178 


4. Effect of different initial separations as a 
substitute for the eccentricity 


others, derived quantities like v^ve and are relatively 
robust with respect to the unphysical eccentricity. 


Ejecta properties computed in our simulations devi¬ 
ate from those of hypothetical genuinely circular mergers 
due to unphysical eccentricities inherent in initial data. 
Specifically, different orbital and approaching velocities 
at the tidal disruption lead to deviations of characteristic 
quantities on the order of the eccentricity. Although this 
error can be eliminated by iterative eccentricity reduction 
[921 1146j , it is demanding to reduce the eccentricities for 
all the models considered in this study. 


To estimate systematic errors associated with unphysi¬ 
cal eccentricities, we instead compare results obtained by 
models with different values of mo^o with e ~ 0.01-0.02 
for APR4-Q3a75. These models should merge at a dif¬ 
ferent true anomaly (angle measured from the periapsis), 
admitting that it is very difficult to quantify this state¬ 
ment in numerical simulations. Therefore, the results 
will give us an idea of errors associated with eccentrici¬ 
ties. As shown in Table VIII ejecta quantities like Mej, 
Tej, and Pej fluctuate within ±2.5% as expected from the 
value of the eccentricity. Because the increase/decrease 
of a single ejecta quantity is accompanied by that of the 


5. Comment on seemingly spurious high-velocity 
ejecta 

In Sec. [IIIB 1[ a small amount of unbound material 
is found to be ejected with a large velocity. We re¬ 
gard this component as an artifact, because the amount 
of high-velocity ejecta does not converge even approxi¬ 
mately with respect to grid resolutions, admitting that 
ejecta cannot be decomposed unambiguously into physi¬ 
cal and unphysical components unless reliable extrapola¬ 
tion to the continuum limit is performed. We speculate 
that this high-velocity ejecta is created by the artificial 
atmosphere and finite grid resolutions. They induce un¬ 
physical shocks at the stellar surface during the inspiral 
phase, and tenuous material continuously flows out from 
the inner edge of the neutron star. This artificial out¬ 
flow is accumulated around the black hole and forms a 
small unphysical disk during the inspiral phase, whereas 
some of this disk may be supplied by the neutron star 
in a physical manner after the onset of mass shedding. 
A fraction of this unphysical disk is ejected impulsively 
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TABLE VII. The same as Table m but for runs with different values of Tth for selected models. The fiducial value of Tth 
adopted in the body text is 1.8. All the simulations are performed with N = 60. 


Lth 

'^r>rAH [-^©] 

Afbd[Al0] 

fofejlMo] 

Tej (erg) 

P.AMq] 

“^ave 

rej 

APR4-Q3a75 

1.6 

0.196 

0.184 

0.0116 

5.95 X 10™ 

2.13 X 10"® 

0.240 

0.184 

1.8 

0.194 

0.182 

0.0111 

5.48 X 10®° 

2.08 X 10“® 

0.235 

0.187 

2.0 

0.192 

0.184 

0.0082 

3.69 X 10®° 

1.51 X 10“® 

0.225 

0.184 

H4-Q3a75 

1.6 

0.326 

0.280 

0.0455 

2.41 X 10®" 

9.34 X 10"® 

0.243 

0.205 

1.8 

0.331 

0.285 

0.0454 

2.30 X 10®" 

9.22 X 10“® 

0.238 

0.203 

2.0 

0.334 

0.287 

0.0464 

2.35 X 10®" 

9.28 X 10“® 

0.238 

0.200 

APR4-Q5a75 

1.6 

0.064 

0.053 

0.0108 

6.75 X 10®° 

8.82 X lO""" 

0.265 

0.082 

1.8 

0.068 

0.059 

0.0084 

4.77 X 10®° 

8.30 X 10“^ 

0.252 

0.099 

2.0 

0.072 

0.063 

0.0089 

5.08 X 10®° 

9.80 X lO”"" 

0.253 

0.110 

H4-Q5a75 

1.6 

0.314 

0.261 

0.0529 

3.74 X 10®" 

1.14 X 10"® 

0.281 

0.216 

1.8 

0.316 

0.266 

0.0502 

3.33 X 10®" 

1.10 X 10“^ 

0.272 

0.220 

2.0 

0.320 

0.269 

0.0516 

3.34 X 10®" 

1.13 X 10“^ 

0.269 

0.220 


TABLE VIII. The same as Table [nl| but with a different initial angular velocity, mollo, for APR4-Q3a75. The fiducial value of 
mofl adopted in the body text is 0.036. All the simulations are performed with N = 60. 


mollo 

'^r>rAH [-^o] 

AfbdjAf©] 


Tej (erg) 

P^AMq] 

"^ave 

dej 

0.036 

0.194 

0.182 

0.0111 

5.48 X 10®° 

2.08 X 10"® 

0.235 

0.187 

0.034 

0.206 

0.194 

0.0112 

5.59 X 10®° 

2.11 X 10"° 

0.236 

0.188 

0.032 

0.203 

0.192 

0.0109 

5.38 X 10®° 

2.02 X 10"° 

0.235 

0.186 

0.030 

0.205 

0.194 

0.0114 

5.64 X 10®° 

2.12 X 10"° 

0.236 

0.186 


during the merger (due possibly to the tidal torque ex¬ 
erted by the neutron star) with a large velocity reflecting 
the large escape velocity of black holes, i.e., the speed of 
light. 


Appendix B: Analytic estimate of ejecta opening 
angle 


As we discussed in Sec. |IIIB 2[ the ejecta geometry 
may be characterized by the opening angle in the merid¬ 
ional plane, and that in the meridional plane, Ogj. 
Before looking at numerical results, it is instructive to 
estimate these angles by analytic arguments for compar¬ 
isons. These estimates help us to distinguish between 
expected and unexpected features. 

Allowing more than one revolution, the opening angle 
of the dynamical ejecta in the equatorial plane should be 
given by 


where p is the average stellar rest-mass density, which is 
determined by the equation of state. On the other hand, 
Ptd should be given by 



(B3) 


where Eq. (36) is used and spin-induced corrections are 


temporarily neglected. This suggests that the depen¬ 
dence of i^ej on the equation of state is weak, because p 
cancels. This expression also suggests that (pei is smaller 
for a larger mass ratio, but the expected change is less 
than 10% between Q = 3 and 7. Prograde black-hole 
spins will decrease tpej, because the orbital frequency 
around a Kerr black hole is given by |SH] 


Ok = 




BH 


r3/2 




(B4) 


Pej 


2tt 


ttd 


td 


(Bl) 


where ttd is the time scale of tidal disruption and Ptd is 
the orbital period at the tidal disruption radius, rtd (see 
Eq. (36)). On one hand, ttd may be given approximately 
by the sound crossing time tgc of the neutron star as 


1 

ttd ~ tgc OC — 




(B2) 


and thus Ptd increases as y increases. 

The opening angle of the dynamical ejecta in the 
meridional plane, Pgj, is determined by the ratio of the 
velocity perpendicular to the orbital plane, 'Uj_, to that 
in the equatorial direction, z)||, as Pgj « arctan(i;j_/z)||) « 
Vd_/v\\. This value should be given by the ratio of the 
neutron-star radius perpendicular to the orbital plane to 
the tidal disruption radius, rtd [see Eq. (36)]. Thus, the 
dependence of 0ej on the equation of state will be weak 
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again, because both uy and z;_l should scale linearly with 
-Rns- Dependence on the mass ratio is expected to be 


0 ej oc inherited from rtd, but the expected change 

is only 25% between Q = 3 and 7. The spin will not 
modify the value of 0ej • 
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